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
>