Estimating the Binomial Probit Model
EC 313  Spring 2000

SetDirectory["Keewaydin:Rumelihisari:Stata:Stata Data"];

Data from Stata's auto.dta:74 observations on automobile models.
y==1 if car is foreign,==0 if domestic
X variables are {price,mpg,weight} of each car

y=ReadList["y.dat"];
x=ReadList["x.dat",{Number,Number,Number}];

<< Statistics`MultiDescriptiveStatistics`

Mean[y] // N

30% of the cars in the sample are foreign makes

Mean[x] // N

The mean price (1987 dollars), mpg, and weight (lbs)

StandardDeviation[x] // N

CorrelationMatrix[x] // N // MatrixForm

Price and MPG are negatively correlated; price and weight are positively correlated, 
and not surprisingly there is a very large negative correlation between MPG and weight.

The Binomial Probit Model

If we are trying to determine whether a particular characteristic is present or absent, whether a particular event  happened or not (or is likely to happen), or estimate the probability that it will happen, we must be careful in 'predicting' the event (or its probability) from other measurable characteristics. Usually, linear regression would be used to 'fit' an outcome variable to a set of explanatory variables. But in this case, linear regression is not appropriate,s ince the left hand side variable is a probability, confined to the unit interval. A regression equation will predict a less-than-zero or greater-than-unity probability for extreme values of the explanatory variables. 
      An alternative formulation of this estimation problem  makes use of the Probit model, in which an index value (akin to the predicted value from least squares regression) is fed through the cumulative distribution function of the Normal distribution. No matter how large or small a value is chosen, the CDF of that value must range between zero (for a completely unlikely event) and one (for an almost certain event). Thus, the probit model can express a binary, or binomial outcome in terms of a nonlinear function of a set of explanatory variables. The estimation of this model is more complicated than regression, since there is no closed form solution for the parameters. But it can be done by maximizing the likelihood function of the sample data, where that function takes the Normal CDF into account. 
      The Probit function defined below takes a vector y, which should contain only (0,1) values,  and a matrix X of explanatory variables of any size. The number of rows of the y vector must equal the number of rows of the X matrix. A constant term is provided, and a "linear probability model" (regression of y on augmented X) is used to generate starting values for the nonlinear optimization routine.
      The LLF function defines the scalar likelihood of the sample, which is the sum of logarithms of values generated from each observation for the current set of parameters. Since the built-in numerical optimization routine is FindMinimum, minus the likelihood is defined so that a minimum of this function will represent the maximized value of the actual likelihood function. 
      After the optimum is located, the estimated coefficients are used to generate predicted values, defined as the CDF of the index value generated from X and the estimated coefficients. Unlike the original data, which take on only 0 and 1 values, the predicted values will generally lie withn the unit interval.

BeginPackage["ProbitPack`"];
Probit::usage = "Probit[y,x] estimates a bivariate probit model by \
maximum likelihood."
LogLikelihood::usage = "LogLikelihood is the maximized value of log \
likelihood."
ProbitCoefficients::usage = "ProbitCoefficients are the point \
estimates." 
yhat::usage = "yhat are the predicted probabilities of success for \
each observation."
Needs["Statistics`LinearRegression`"];
Needs["Statistics`NormalDistribution`"];
Needs["Statistics`Common`RegressionCommon`"];
Begin["`Private`"];  

\(Probit[y_List, x_List] := 
      Block[{xaug, iota, z, betahat, \ loglike, bhl, startval, probitarg, 
          ProbitOutput, 
estbeta, theta, 
          ndist}, 
n = \(Dimensions[
              y]\)[\([1]\)]; 
 (*\ 
          get\ starting\ values\ from\ linear\ probability\ model\ *) 
 (*\ create\ a\ constant\ term\ *) 
	
        iota = Table[1, {n}]; \ 
 (*\ 
          augment\ the\ X\ matrix\ with\ the\ constant\ term\ *) 
	
        xaug = Transpose[Join[{iota}, Transpose[x]]]; 
 (*\ 
          regress\ y\ on\ Xaugmented\ *) 
	
        z = DesignedRegress[xaug, y]; 
 (*\ 
          extract\ the\ betahat\ point\ estimates\ *) 
	
        betahat = \(Transpose[\((ParameterTable /. 
                    z)\)[\([1]\)]]\)[\([1]\)]\ ; 
 (*\ 
          generate\ a\ list\ of\ parameters\ *) 
	
        bhl = Table[
            Global`b[i], {i, Length[betahat]}]; 
 (*\ 
          
          thread\ the\ list\ of\ starting\ values\ through\ the\ parameter\ list\ *) 
	
        startval = MapThread[List, {bhl, betahat}]; 
 (*\ 
          use\ FindMinimum\ to\ locate\ minus\ the\ min\ likelihood\ *) 
	
        probitarg = 
          Append[Prepend[startval, 
              Hold[LLF[y, xaug, bhl, n]]], {Method -> 
                Newton}]; 
	
        ProbitOutput = 
          Apply[FindMinimum, ReleaseHold[probitarg]]; 
	
        LogLikelihood = \(-ProbitOutput[\([1]\)]\); 
	
        estbeta = bhl /. ProbitOutput[\([2]\)]; 
	
        theta = xaug . estbeta; 
	
        ndist = NormalDistribution[0, 1]; 
	
        yhat = Table[0, {n}]; 
	
        Do[yhat[\([i]\)] = CDF[ndist, theta[\([i]\)]], {i, 1, 
            n}]; 
	
        ProbitCoefficients = 
          ProbitOutput[\([2]\)] // 
            TableForm
]\ 
\ \)
LLF[y_, x_, beta_, n_] :=  
      Block[{theta, ndist, i,  elelf}, 
theta = 
          x . beta; 
elelf = 
          0.0; 
ndist = 
          NormalDistribution[0, 1]; 
Do[
          If[y[[i]] == 1, 
elelf -= 
              Log[CDF[ndist, theta[[i]]]], 
            elelf  -= Log[CDF[ndist, -theta[[i]]]]], {i, 1, 
            n}]; 
elelf]

End[]; 
EndPackage[];

Probit[y, x]

This indicates that price, b[2], has a very small positive effect on the probability that a car is foreign, while both higher MPG and higher weight make it more likely that a car is domestic. The former result may be surprising, since you would expect foreign cars to get better mileage. They do,but these are partial derivatives, so that holding price and weight constant, the ∂Pr[foreign]/ ∂MPG is negative. Note that the estimated coefficients are not, in themselves, estimates of how much the probability will change for a one unit change in the explanatory variable; they are estimates of how much the unobserved 'index variable' will change. The chain rule must then be used to evaluate the effect on the probability, evaluated at a particular point in the X space. Thus, unlike regression, where the effects of an explanatory variable are constant, they vary continuously in the probit model.

LogLikelihood

Mean[yhat]

ListPlot[yhat]

dompred = Extract[yhat, Position[y, 0]];

forpred = Extract[yhat, Position[y, 1]];

Mean[dompred]

Median[dompred]

Mean[forpred]

Median[forpred]

Note that most domestic cars are correctly classified as such by the model; it is much less successful with foreign cars.

The probit model is commonly applied in decision making processes, such as extension of credit or hiring a worker, where a number of disparate factors may be observed and a 'yes/no' decision must be reached. It is also widely employed in applied econometric work to evaluate the factors that are relevant in consumers' choices--e.g. in evaluating the success of a marketing campaign, or voting patterns.

Acknowledgement: I am very grateful to David Belsley for his expert assistance with Mathematica syntax.