MathGroup Archive 1998

[Date Index] [Thread Index] [Author Index]

Search the Archive

FindRoot, a big Jacobian, and bivariate probit model

  • To: mathgroup at smc.vnet.net
  • Subject: [mg14107] FindRoot, a big Jacobian, and bivariate probit model
  • From: "ELLIS, Luci (FS)" <EllisL at rba.gov.au>
  • Date: Fri, 25 Sep 1998 03:15:31 -0400
  • Sender: owner-wri-mathgroup at wolfram.com

> Dear MathGroup,
> I am currently trying to write Mathematica code for estimating a bivariate
> probit model by Maximum Likelihood and have run into some problems (to say
> the least).  I am not getting any error messages, but even on a toy
> problem
> (3 explanatory variables, 20 data points), Mathematica fails to give an
> answer in any reasonable amount of time.  
> 
> (For those familiar with FindRoot but not probit models, the short
> explanation is that they are models where the dependent variable takes the
> values 0 or 1, that is,  y = 1 if b.X + e > 0 and 0 otherwise, where e is
> a
> standard normally-distributed disturbance.  In the bivariate case, you
> have
> two y-variables, and the errors for each of them may be correlated, ie
> distributed bivariate normal)
> 
> I have analytic expressions for the first-order conditions and the
> Jacobian
> (second derivatives of the log-likelihood, in that case).  However, the
> Jacobian takes a L-O-N-G time to evaluate because it involves the CDF of a
> bivariate normal distribution, and the parameters are not known, which
> means
> that they have to be substituted into this enormous expression.  I have
> already tried NOT giving Mathematica the Jacobian explicitly, but then it
> asks me for one because it can't get the analytic derivatives of these
> expressions.  How can I speed this up so that the code is usable?
> 
> Thanks in advance,
> Luci Ellis
> 
> The function currently looks like this: (note I haven't Modulised all the
> local variables yet, so I can look at them.  You'll have to copy and paste
> this into Mathematica to see it properly.)
> (****start from next line ****)
> Cell[BoxData[
>     RowBox[{
>       RowBox[{
>         RowBox[{"EstimateBivariateProbitModel", "[", 
>           RowBox[{
>             RowBox[{"y1", ":", 
>               RowBox[{"{", 
>                 RowBox[{
>                   RowBox[{"(", 
>                     RowBox[{"0", "|", "1"}], ")"}], ".."}], "}"}]}], ",", 
>             RowBox[{"y2", ":", 
>               RowBox[{"{", 
>                 RowBox[{
>                   RowBox[{"(", 
>                     RowBox[{"0", "|", "1"}], ")"}], ".."}], "}"}]}], ",", 
>             RowBox[{"x_", "?", "MatrixQ"}], ",", "opts___Rule"}], "]"}], 
>         "/;", 
>         RowBox[{"(", 
>           RowBox[{
>             RowBox[{"(", 
>               RowBox[{
>                 RowBox[{"Length", "[", "y1", "]"}], "==", 
>                 RowBox[{"Length", "[", "y2", "]"}], "==", 
>                 RowBox[{"Length", "[", "x", "]"}]}], ")"}], " ", "&&", "
> ", 
>             RowBox[{"And", " ", "@@", 
>               RowBox[{"(", " ", 
>                 RowBox[{"NumericQ", "/@", 
>                   RowBox[{"Flatten", "[", "x", "]"}]}], ")"}]}]}],
> ")"}]}], 
>       ":=", "\n", "\t\t", 
>       RowBox[{"Module", "[", 
>         RowBox[{
>           RowBox[{"{", "}"}], ",", "\n", "\t", 
>           RowBox[{"With", "[", 
>             RowBox[{
>               RowBox[{"{", 
>                 RowBox[{
>                   RowBox[{"\[CapitalBeta]1", "=", 
>                     RowBox[{"ToExpression", "[", 
>                       RowBox[{
>                         RowBox[{
>                           RowBox[{"\"\<\[Beta]\>\"", "<>", "#"}], " ",
> "&"}], 
>                         " ", "/@", " ", 
>                         RowBox[{"(", 
>                           RowBox[{"ToString", " ", "/@", " ", 
>                             RowBox[{"Range", "[", 
>                               RowBox[{
>                                 RowBox[{"Dimensions", "[", "x", "]"}], 
>                                 "\[LeftDoubleBracket]", "2", 
>                                 "\[RightDoubleBracket]"}], "]"}]}],
> ")"}]}],
> 
>                       "]"}]}], ",", 
>                   RowBox[{"\[CapitalBeta]2", "=", 
>                     RowBox[{"ToExpression", "[", 
>                       RowBox[{
>                         RowBox[{
>                           RowBox[{"\"\<\[Gamma]\>\"", "<>", "#"}], " ", 
>                           "&"}], " ", "/@", " ", 
>                         RowBox[{"(", 
>                           RowBox[{"ToString", " ", "/@", " ", 
>                             RowBox[{"Range", "[", 
>                               RowBox[{
>                                 RowBox[{"Dimensions", "[", "x", "]"}], 
>                                 "\[LeftDoubleBracket]", "2", 
>                                 "\[RightDoubleBracket]"}], "]"}]}],
> ")"}]}],
> 
>                       "]"}]}], ",", " ", 
>                   RowBox[{"q1", "=", 
>                     RowBox[{
>                       RowBox[{"2", " ", "y1"}], "-", "1"}]}], ",", 
>                   RowBox[{"q2", " ", "=", 
>                     RowBox[{
>                       RowBox[{"2", " ", "y2"}], "-", "1"}]}]}], "}"}],
> ",", 
>               "\n", "\t\t", 
>               RowBox[{
>                 RowBox[{"z1", "=", 
>                   RowBox[{"x", ".", "\[CapitalBeta]1"}]}], ";", "\t", 
>                 RowBox[{"z2", "=", 
>                   RowBox[{"x", ".", "\[CapitalBeta]2"}]}], ";", 
>                 RowBox[{"\[Rho]star", "=", " ", 
>                   RowBox[{"q1", "*", "q2", "*", "\[Rho]"}]}], ";", 
>                 RowBox[{"w1", "=", 
>                   RowBox[{"q1", "*", "z1"}]}], ";", " ", 
>                 RowBox[{"w2", "=", 
>                   RowBox[{"q2", "*", "z2"}]}], ";", "\n", "\t\t\t", 
>                 RowBox[{"g1", "=", 
>                   RowBox[{"Evaluate", "[", 
>                     RowBox[{
>                       RowBox[{"\[Phi]", "[", "w1", "]"}], "*", 
>                       RowBox[{"\[CapitalPhi]", "[", 
>                         FractionBox[
>                           RowBox[{"w2", "-", 
>                             RowBox[{"\[Rho]star", " ", "w1"}]}], 
>                           SqrtBox[
>                             RowBox[{"(", 
>                               RowBox[{"1", "-", 
>                                 SuperscriptBox["\[Rho]star", "2"]}],
> ")"}]]], 
>                         "]"}]}], "]"}]}], ";", "\n", "\t\t\t", 
>                 RowBox[{"g2", "=", 
>                   RowBox[{"Evaluate", "[", 
>                     RowBox[{
>                       RowBox[{"\[Phi]", "[", "w2", "]"}], "*", 
>                       RowBox[{"\[CapitalPhi]", "[", 
>                         FractionBox[
>                           RowBox[{"w1", "-", 
>                             RowBox[{"\[Rho]star", " ", "w2"}]}], 
>                           SqrtBox[
>                             RowBox[{"(", 
>                               RowBox[{"1", "-", 
>                                 SuperscriptBox["\[Rho]star", "2"]}],
> ")"}]]], 
>                         "]"}]}], "]"}]}], ";", "\t\t\t", "\n", "\t\t\t", 
>                 RowBox[{"triple", "=", 
>                   RowBox[{"Transpose", "[", 
>                     RowBox[{"{", 
>                       RowBox[{"w1", ",", "w2", ",", "\[Rho]star"}], "}"}],
> 
>                     "]"}]}], ";", "\n", "\t\t\t", 
>                 RowBox[{"\[Phi]2", "=", 
>                   RowBox[{
>                     RowBox[{
>                       RowBox[{"PDF", "[", 
>                         RowBox[{
>                           RowBox[{"MultinormalDistribution", "[", 
>                             RowBox[{
>                               RowBox[{"{", 
>                                 RowBox[{"0", ",", "0"}], "}"}], ",", 
>                               TagBox[
>                                 RowBox[{"{", 
>                                   RowBox[{
>                                     RowBox[{"{", 
>                                       RowBox[{"1", ",", 
>                                         RowBox[{
>                                         StyleBox["#",
>                                         StyleBoxAutoDelete->True,
>                                         PrintPrecision->1], 
>                                         "\[LeftDoubleBracket]", "3", 
>                                         "\[RightDoubleBracket]"}]}],
> "}"}], 
>                                     ",", 
>                                     RowBox[{"{", 
>                                       RowBox[{
>                                         RowBox[{
>                                         StyleBox["#",
>                                         StyleBoxAutoDelete->True,
>                                         PrintPrecision->1], 
>                                         "\[LeftDoubleBracket]", "3", 
>                                         "\[RightDoubleBracket]"}], ",", 
>                                         "1"}], "}"}]}], "}"}],
>                                 Short]}], "]"}], ",", 
>                           RowBox[{"{", 
>                             RowBox[{
>                               RowBox[{
>                               "#", "\[LeftDoubleBracket]", "1", 
>                                 "\[RightDoubleBracket]"}], ",", 
>                               RowBox[{
>                               "#", "\[LeftDoubleBracket]", "2", 
>                                 "\[RightDoubleBracket]"}]}], "}"}]}],
> "]"}],
> 
>                       "&"}], " ", "/@", " ", "triple"}]}], ";", "\n", 
>                 "\t\t\t", 
>                 RowBox[{"\[CapitalPhi]2", "=", 
>                   RowBox[{
>                     RowBox[{
>                       RowBox[{"CDF", "[", 
>                         RowBox[{
>                           RowBox[{"MultinormalDistribution", "[", 
>                             RowBox[{
>                               RowBox[{"{", 
>                                 RowBox[{"0", ",", "0"}], "}"}], ",", 
>                               TagBox[
>                                 RowBox[{"{", 
>                                   RowBox[{
>                                     RowBox[{"{", 
>                                       RowBox[{"1", ",", 
>                                         RowBox[{
>                                         StyleBox["#",
>                                         StyleBoxAutoDelete->True,
>                                         PrintPrecision->1], 
>                                         "\[LeftDoubleBracket]", "3", 
>                                         "\[RightDoubleBracket]"}]}],
> "}"}], 
>                                     ",", 
>                                     RowBox[{"{", 
>                                       RowBox[{
>                                         RowBox[{
>                                         StyleBox["#",
>                                         StyleBoxAutoDelete->True,
>                                         PrintPrecision->1], 
>                                         "\[LeftDoubleBracket]", "3", 
>                                         "\[RightDoubleBracket]"}], ",", 
>                                         "1"}], "}"}]}], "}"}],
>                                 Short]}], "]"}], ",", 
>                           RowBox[{"{", 
>                             RowBox[{
>                               RowBox[{
>                               "#", "\[LeftDoubleBracket]", "1", 
>                                 "\[RightDoubleBracket]"}], ",", 
>                               RowBox[{
>                               "#", "\[LeftDoubleBracket]", "2", 
>                                 "\[RightDoubleBracket]"}]}], "}"}]}],
> "]"}],
> 
>                       "&"}], " ", "/@", " ", "triple"}]}], ";", "\n", 
>                 "\t\t\t", 
>                 RowBox[{"\[CapitalPhi]2sq", "=", 
>                   RowBox[{"\[CapitalPhi]2", "^", "2"}]}], ";", "\n", 
>                 "\t\t\t", 
>                 RowBox[{"\[Phi]\[CapitalPhi]", "=", 
>                   RowBox[{"\[Phi]2", "/", "\[CapitalPhi]2"}]}], ";", "\n",
> 
>                 "\t\t\t", 
>                 RowBox[{"\[Delta]", "=", 
>                   RowBox[{"1", "/", 
>                     SqrtBox[
>                       RowBox[{"(", 
>                         RowBox[{"1", "-", 
>                           RowBox[{"\[Rho]star", "^", "2"}]}], ")"}]]}]}], 
>                 ";", "\n", "\t\t\t", 
>                 RowBox[{"v1", "=", 
>                   RowBox[{"\[Delta]", 
>                     RowBox[{"(", 
>                       RowBox[{"w2", "-", 
>                         RowBox[{"\[Rho]star", " ", "w1"}]}], ")"}]}]}],
> ";",
> 
>                 "\n", "\t\t\t", 
>                 RowBox[{"v2", "=", 
>                   RowBox[{"\[Delta]", " ", 
>                     RowBox[{"(", 
>                       RowBox[{"w1", "-", 
>                         RowBox[{"\[Rho]star", " ", "w2"}]}], ")"}]}]}],
> ";",
> 
>                 "\n", "\t\t\t", 
>                 RowBox[{"focs", "=", 
>                   RowBox[{"Evaluate", "[", 
>                     RowBox[{"Thread", "[", 
>                       RowBox[{
>                         RowBox[{"Join", "[", 
>                           RowBox[{
>                             RowBox[{
>                               RowBox[{"(", 
>                                 RowBox[{"q1", " ", 
>                                   RowBox[{"g1", "/", "\[CapitalPhi]2"}]}],
> 
>                                 ")"}], ".", "x"}], ",", 
>                             RowBox[{
>                               RowBox[{"(", 
>                                 RowBox[{"q2", " ", 
>                                   RowBox[{"g2", "/", "\[CapitalPhi]2"}]}],
> 
>                                 ")"}], ".", "x"}], ",", 
>                             RowBox[{"{", 
>                               RowBox[{
>                                 RowBox[{"(", 
>                                   RowBox[{"q1", " ", 
>                                     RowBox[{"q2", "/",
> "\[CapitalPhi]2"}]}],
> 
>                                   ")"}], ".", "\[Phi]2"}], " ", "}"}]}], 
>                           "]"}], "==", "0"}], "]"}], "]"}]}], ";", "\n", 
>                 "\t\t\t", 
>                 RowBox[{"\[ScriptD]b1b1", "=", 
>                   RowBox[{
>                     RowBox[{"Transpose", "[", "x", "]"}], ".", 
>                     RowBox[{"(", 
>                       RowBox[{"x", "*", 
>                         RowBox[{"(", " ", 
>                           RowBox[{
>                             RowBox[{"(", 
>                               RowBox[{
>                                 RowBox[{"-", "w1"}], " ", 
>                                 RowBox[{"g1", "/", "\[CapitalPhi]2"}]}], 
>                               ")"}], "-", 
>                             RowBox[{"(", 
>                               RowBox[{
>                               "\[Rho]star", "  ", "\[Phi]\[CapitalPhi]"}],
> 
>                               ")"}], " ", "-", 
>                             RowBox[{"(", 
>                               RowBox[{
>                                 RowBox[{"g1", "^", "2"}], "/", 
>                                 "\[CapitalPhi]2sq"}], ")"}]}], ")"}]}], 
>                       ")"}]}]}], ";", "\n", "\t\t\t", 
>                 RowBox[{"\[ScriptD]b1b2", "=", 
>                   RowBox[{
>                     RowBox[{"Transpose", "[", 
>                       RowBox[{"q1", " ", "q2", " ", "x"}], "]"}], ".", 
>                     RowBox[{"(", 
>                       RowBox[{"x", "*", " ", 
>                         RowBox[{"(", 
>                           RowBox[{"\[Phi]\[CapitalPhi]", "-", 
>                             RowBox[{"(", 
>                               RowBox[{"g1", " ", 
>                                 RowBox[{
>                                 "g2", " ", "/", "\[CapitalPhi]2sq"}]}], 
>                               ")"}]}], ")"}]}], ")"}]}]}], ";", "\n", 
>                 "\t\t\t", 
>                 RowBox[{"\[ScriptD]b2b2", "=", 
>                   RowBox[{
>                     RowBox[{"Transpose", "[", "x", "]"}], ".", 
>                     RowBox[{"(", 
>                       RowBox[{"x", "*", 
>                         RowBox[{"(", " ", 
>                           RowBox[{
>                             RowBox[{"(", 
>                               RowBox[{
>                                 RowBox[{"-", "w1"}], " ", 
>                                 RowBox[{"g1", "/", "\[CapitalPhi]2"}]}], 
>                               ")"}], "-", 
>                             RowBox[{"(", 
>                               RowBox[{
>                               "\[Rho]star", "  ", "\[Phi]\[CapitalPhi]"}],
> 
>                               ")"}], " ", "-", 
>                             RowBox[{"(", 
>                               RowBox[{
>                                 RowBox[{"g1", "^", "2"}], "/", 
>                                 "\[CapitalPhi]2sq"}], ")"}]}], ")"}]}], 
>                       ")"}]}]}], ";", "\n", "\t\t\t", 
>                 RowBox[{"\[ScriptD]b1b\[Rho]", "=", 
>                   RowBox[{"{", 
>                     RowBox[{
>                       RowBox[{"Transpose", "[", 
>                         RowBox[{"q2", "  ", "x"}], "]"}], ".", 
>                       RowBox[{"(", 
>                         RowBox[{"\[Phi]\[CapitalPhi]", 
>                           RowBox[{"(", 
>                             RowBox[{
>                               RowBox[{
>                               "\[Rho]star", " ", "\[Delta]", " ", "v1"}], 
>                               "-", "w1", " ", "-", 
>                               RowBox[{"g1", "/", "\[CapitalPhi]2"}]}], 
>                             ")"}]}], ")"}]}], "}"}]}], ";", "\n",
> "\t\t\t", 
>                 RowBox[{"\[ScriptD]b2b\[Rho]", "=", 
>                   RowBox[{"{", 
>                     RowBox[{
>                       RowBox[{"Transpose", "[", 
>                         RowBox[{"q1", "  ", "x"}], "]"}], ".", 
>                       RowBox[{"(", 
>                         RowBox[{"\[Phi]\[CapitalPhi]", 
>                           RowBox[{"(", 
>                             RowBox[{
>                               RowBox[{
>                               "\[Rho]star", " ", "\[Delta]", " ", "v2"}], 
>                               "-", "w2", "-", 
>                               RowBox[{"g2", "/", "\[CapitalPhi]2"}]}], 
>                             ")"}]}], ")"}]}], "}"}]}], ";", "\n",
> "\t\t\t", 
>                 RowBox[{"\[ScriptD]\[Rho]\[Rho]", "=", 
>                   RowBox[{"{", 
>                     RowBox[{"{", 
>                       RowBox[{"\[Phi]\[CapitalPhi]", ".", 
>                         RowBox[{"(", 
>                           RowBox[{
>                             RowBox[{
>                               RowBox[{"\[Delta]", "^", "2"}], "  ", 
>                               "\[Rho]star", " ", 
>                               RowBox[{"(", 
>                                 RowBox[{"1", "-", 
>                                   RowBox[{
>                                     RowBox[{"\[Delta]", "^", "2"}], " ", 
>                                     RowBox[{"(", 
>                                       RowBox[{
>                                         RowBox[{"w1", "^", "2"}], "+", 
>                                         RowBox[{"w2", "^", "2"}], "-", 
>                                         RowBox[{
>                                         "2", " ", "\[Rho]star", " ", "w1",
> 
>                                         " ", "w2"}]}], ")"}]}]}], ")"}]}],
> 
>                             "+", 
>                             RowBox[{"\[Delta]", " ", "w1", " ", "w2"}], "
> ",
> 
>                             "-", "\[Phi]\[CapitalPhi]"}], ")"}]}], "}"}], 
>                     "}"}]}], ";", "\n", "\t\t\t", 
>                 RowBox[{"jac", "=", 
>                   RowBox[{"BlockMatrix", "[", 
>                     RowBox[{"{", 
>                       RowBox[{
>                         RowBox[{"{", 
>                           RowBox[{
>                           "\[ScriptD]b1b1", ",", "\[ScriptD]b1b2", ",", 
>                             RowBox[{
>                             "Transpose", "[", "\[ScriptD]b1b\[Rho]",
> "]"}]}], 
>                           "}"}], ",", 
>                         RowBox[{"{", 
>                           RowBox[{
>                             RowBox[{
>                             "Transpose", "[", "\[ScriptD]b1b2", "]"}],
> ",", 
>                             "\[ScriptD]b2b2", ",", 
>                             RowBox[{
>                             "Transpose", "[", "\[ScriptD]b2b\[Rho]",
> "]"}]}], 
>                           "}"}], ",", 
>                         RowBox[{"{", 
>                           RowBox[{
>                           "\[ScriptD]b1b\[Rho]", ",",
> "\[ScriptD]b2b\[Rho]",
> 
>                             ",", "\[ScriptD]\[Rho]\[Rho]"}], "}"}]}],
> "}"}],
> 
>                     "]"}]}], ";", "\n", "\t", 
>                 RowBox[{"res", "=", "\t", 
>                   RowBox[{"FindRoot", "[", 
>                     RowBox[{
>                       RowBox[{"Evaluate", "[", "focs", "]"}], ",", 
>                       RowBox[{"Evaluate", "[", 
>                         RowBox[{"Sequence", " ", "@@", 
>                           RowBox[{"Transpose", "[", 
>                             RowBox[{"{", 
>                               RowBox[{
>                                 RowBox[{"Join", "[", 
>                                   RowBox[{
>                                   "\[CapitalBeta]1", ",",
> "\[CapitalBeta]2",
> 
>                                     ",", 
>                                     RowBox[{"{", "\[Rho]", "}"}]}], "]"}],
> 
>                                 ",", 
>                                 RowBox[{"Table", "[", 
>                                   RowBox[{"0.01", ",", 
>                                     RowBox[{"{", 
>                                       RowBox[{
>                                         RowBox[{
>                                         "Length", "[", "\[CapitalBeta]1", 
>                                         "]"}], "+", 
>                                         RowBox[{
>                                         "Length", "[", "\[CapitalBeta]2", 
>                                         "]"}], "+", "1"}], "}"}]}],
> "]"}]}],
> 
>                               "}"}], "]"}]}], "]"}], ",", 
>                       RowBox[{"Jacobian", "\[Rule]", "jac"}], ",", 
>                       RowBox[{"Evaluate", "[", 
>                         RowBox[{"FilterOptions", "[", 
>                           RowBox[{"FindRoot", ",", "opts", " ", ",", 
>                             RowBox[{"Sequence", " ", "@@", 
>                               RowBox[{"Options", "[", "FindRoot",
> "]"}]}]}],
> 
>                           " ", "]"}], "]"}]}], " ", "]"}]}], ";", " ",
> "\n",
> 
>                 "\t\t\t\t", 
>                 RowBox[{"{", "res", "}"}]}]}], "\n", "\t", "]"}]}],
> "]"}]}]], 
>   "Input"]
> 
> ________________________________________________________
> Luci Ellis
> x8786
> Payments Policy Department
> ellisl at rba.gov.au
> 


  • Prev by Date: finite element analysis in mathematica
  • Next by Date: Re: Fourier transform
  • Previous by thread: finite element analysis in mathematica
  • Next by thread: Re: How to import DTA files conatining comment lines starting with '#'