(* ::Package:: *) (* : Title : Logit and Probit Models*) (* : Author :L. Aylin Aktukun*) (* : Mathematica Version : 6.0*) (* : Date : January, 2008*) (* : Summary : This package is designed to calculate the maximum likelihood estimates of parameters, the asymptotic standard errors,t-stats and p-values for logit and probit models.*) BeginPackage["BinaryResponse`"] Logit::usage="calculates the maximum likelihood estimates of parameters, the asymptotic standard errors,t-stats and p-values for logit model"; Probit::usage="calculates the maximum likelihood estimates of parameters, the asymptotic standard errors,t-stats and p-values for probit model"; Begin["`Private`"] Needs["LinearRegression`"] Off[FindRoot::"precw"] Off[FindRoot::"bddir"] Off[General::"spell1"] Logit[data_]:=Module[{n,k,dep,inds,rr,regress,paramEstim,ones,F,desMatrix,betas,xbeta,fonk,equations,solutions,MLEestim,gammaBetas,gammaMatrix,CovMat,stanErrors,tratio,f,pvalues,para}, n=Length[data]; k=Length[Transpose[data]]; dep=Transpose[data][[k]]; inds=Take[Transpose[data],k-1]; rr=Table[Subscript[u,i],{i,1,k-1}]; regress=Regress[data,Flatten[Insert[{1},{rr},2]],rr,RegressionReport->{ParameterTable,BestFitParameters}]; paramEstim=BestFitParameters/.regress; ones=Table[1,{i,1,n}]; F=LogisticDistribution[0,1]; desMatrix=Transpose[Prepend[inds,ones]]; betas=Table[Subscript[b,i],{i,0,k-1}]; xbeta=betas.Transpose[desMatrix]; fonk=dep.Log[CDF[F,xbeta]]+(1-dep).Log[1-CDF[F,xbeta]]; equations=Table[D[fonk,betas[[i]]] ==0,{i,1,k}]; solutions=FindRoot[equations,Transpose[{betas,paramEstim}],WorkingPrecision->10]; MLEestim=Table[betas[[i]]/.solutions,{i,1,k}]; gammaBetas=(PDF[F,xbeta])^2/((CDF[F,xbeta])*(1-CDF[F,xbeta])); gammaMatrix=DiagonalMatrix[gammaBetas/.Table[Subscript[b,i-1]->MLEestim[[i]],{i,1,Length[betas]}]]; CovMat=Inverse[(Transpose[desMatrix].gammaMatrix).desMatrix]; stanErrors=Sqrt[Table[CovMat[[i,i]],{i,1,Length[betas]}]]; tratio=MLEestim/stanErrors; f=StudentTDistribution[n-k]; pvalues=Table[1-CDF[f,Abs[tratio][[i]]],{i,1,Length[betas]}]; para=Table[Subscript["a",i],{i,0,k-1}]; Print["*** Maximum Likelihood Estimates of Parameters for Logit Model ***"]; Print[TableForm[{{"Parameters","MLE Estim.","Asymp. SE.","t Stat","P-Values"},{para,SetPrecision[MLEestim,4],SetPrecision[stanErrors,4],SetPrecision[tratio,3],SetPrecision[pvalues,4] }}]]; ]; Probit[data_]:=Module[{n,k,dep,b,u,inds,rr,regress,paramEstim,ones,F,desMatrix,xbeta,betas,fonk,equations,solutions,MLEestim,gammaBetas,gammaMatrix,CovMat,stanErrors,tratio,f,pvalues,index,para,allparas}, n=Length[data]; k=Length[Transpose[data]]; dep=Transpose[data][[k]]; inds=Take[Transpose[data],k-1]; rr=Table[Subscript[u,i],{i,1,k-1}]; regress=Regress[data,Flatten[Insert[{1},{rr},2]],rr,RegressionReport->{ParameterTable,BestFitParameters}]; paramEstim=BestFitParameters/.regress; ones=Table[1,{i,1,n}]; F=NormalDistribution[0,1]; desMatrix=Transpose[Prepend[inds,ones]]; betas=Table[Subscript[b,i],{i,0,k-1}]; xbeta=betas.Transpose[desMatrix]; fonk=dep.Log[CDF[F,xbeta]]+(1-dep).Log[1-CDF[F,xbeta]]; equations=Table[D[fonk,betas[[i]]] ==0,{i,1,k}]; solutions=FindRoot[equations,Transpose[{betas,paramEstim}],WorkingPrecision->10]; MLEestim=Table[betas[[i]]/.solutions,{i,1,k}]; gammaBetas=(PDF[F,xbeta])^2/((CDF[F,xbeta])*(1-CDF[F,xbeta])); gammaMatrix=DiagonalMatrix[gammaBetas/.Table[Subscript[b,i-1]->MLEestim[[i]],{i,1,Length[betas]}]]; CovMat=Inverse[(Transpose[desMatrix].gammaMatrix).desMatrix]; stanErrors=Sqrt[Table[CovMat[[i,i]],{i,1,Length[betas]}]]; tratio=MLEestim/stanErrors; f=StudentTDistribution[n-k]; pvalues=Table[1-CDF[f,Abs[tratio][[i]]],{i,1,Length[betas]}]; para=Table[Subscript["a",i],{i,0,k-1}]; Print["*** Maximum Likelihood Estimates of Parameters for Probit Model ***"]; Print[TableForm[{{"Parameters","MLE Estim.","Asymp. SE.","t Stat","P-Values"},{para,SetPrecision[MLEestim,4],SetPrecision[stanErrors,4],SetPrecision[tratio,3],SetPrecision[pvalues,4] }}]]; ]; End[] EndPackage[]