Speed up CLR critical value function via reparametrization

[mlondschien]

No description provided.

[mlondschien]

here:

· Discovering benchmarks
· Running 1 total benchmarks (1 commits * 1 environments * 1 benchmarks)
[ 0.00%] ·· Benchmarking existing-py_Users_mlondschien_mambaforge_envs_ivmodels_bin_python
[50.00%] ··· Running (tests.Tests.time_conditional_likelihood_ratio_test_numerical_integration--).
[100.00%] ··· tests.Tests.time_conditional_likelihood_ratio_test_numerical_integration                                                                      ok
[100.00%] ··· ====== ============== ================ ============== =======================
              --                                      data                                 
              ------ ----------------------------------------------------------------------
                n     (1, 1, 0, 0)   (100, 1, 4, 0)   (4, 2, 2, 2)   guggenberger12 (k=10) 
              ====== ============== ================ ============== =======================
               1000     455±8μs         65.5±2ms      1.11±0.03ms          21.5±0.3ms      
              ====== ============== ================ ============== =======================

main:

· Discovering benchmarks
· Running 1 total benchmarks (1 commits * 1 environments * 1 benchmarks)
[ 0.00%] ·· Benchmarking existing-py_Users_mlondschien_mambaforge_envs_ivmodels_bin_python
[50.00%] ··· Running (tests.Tests.time_conditional_likelihood_ratio_test_numerical_integration--).
[100.00%] ··· tests.Tests.time_conditional_likelihood_ratio_test_numerical_integration                                                                      ok
[100.00%] ··· ====== ============== ================ ============== =======================
              --                                      data                                 
              ------ ----------------------------------------------------------------------
                n     (1, 1, 0, 0)   (100, 1, 4, 0)   (4, 2, 2, 2)   guggenberger12 (k=10) 
              ====== ============== ================ ============== =======================
               1000     453±10μs       76.8±0.9ms     1.08±0.02ms         16.1±0.09ms      
              ====== ============== ================ ============== =======================

(ivmodels) (base) code/ivmodels/benchmarks [main] $ 

[mlondschien]

reparametrization is enough:

[ 0.00%] ·· Benchmarking existing-py_Users_mlondschien_mambaforge_envs_ivmodels_bin_python
[50.00%] ··· Running (tests.Tests.time_conditional_likelihood_ratio_test_numerical_integration--).
[100.00%] ··· tests.Tests.time_conditional_likelihood_ratio_test_numerical_integration                                                                      ok
[100.00%] ··· ====== ============== ================ ============== =======================
              --                                      data                                 
              ------ ----------------------------------------------------------------------
                n     (1, 1, 0, 0)   (100, 1, 4, 0)   (4, 2, 2, 2)   guggenberger12 (k=10) 
              ====== ============== ================ ============== =======================
               1000     452±9μs         40.7±2ms      1.08±0.01ms           877±10μs       
              ====== ============== ================ ============== =======================

kleibergen19 script:

(ivmodels) ~/code/ivmodels-simulations [guggenberger12] $ python -m kernprof -lvr experiments/kleibergen19_size_collect.py
Wrote profile results to kleibergen19_size_collect.py.lprof
Timer unit: 1e-06 s

Total time: 0.027289 s
File: /Users/mlondschien/code/ivmodels/ivmodels/tests/conditional_likelihood_ratio.py
Function: conditional_likelihood_ratio_critical_value_function at line 10

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    10                                           @line_profiler.profile
    11                                           def conditional_likelihood_ratio_critical_value_function(
    12                                               p, q, s_min, z, method="numerical_integration", tol=1e-6
    13                                           ):
    14                                               """
    15                                               Approximate the critical value function of the conditional likelihood ratio test.
    16                                           
    17                                               Let
    18                                           
    19                                               .. math: \\Gamma(q-p, p, s_\\mathrm{min}) := 1/2 \\left( Q_{q-p} + Q_p - s_\\mathrm{min} + \\sqrt{ (Q_{q-p} + Q_p - s_\\mathrm{min})^2 + 4 Q_{p} s_\\mathrm{min} } \\right),
    20                                           
    21                                               where :math:`Q_p \\sim \\chi^2(p)` and :math:`Q_{q-p} \\sim \\chi^2(q - p)` are
    22                                               independent chi-squared random variables. This function approximates
    23                                           
    24                                               .. math: Q_{q, p}(z) := \\mathbb{P}[ \\Gamma(q-p, p, s_\\mathrm{min}) > z ]
    25                                           
    26                                               up to tolerance ``tol``.
    27                                           
    28                                               If ``method`` is ``"numerical_integration"``, numerically integrates the formulation
    29                                           
    30                                               .. math: Q_{q, p}(z) = \\mathbb{E}_{B \\sim \\mathrm{Beta}((k - m)/2, m/2)}[ F_{\\chi^2(k)}(z / (1 - a B)) ],
    31                                           
    32                                               where :math:`F_{\\chi^2(k)}` is the cumulative distribution function of a
    33                                               :math:`\\chi^2(k)` distribution, and :math:`a = s_{\\min} / (z + s_{\\min})`. This
    34                                               is Equation (27) of :cite:p:`hillier2009conditional` or Equation (40) of
    35                                               :cite:p:`hillier2009exact`.
    36                                           
    37                                               If ``method`` is ``"power_series"``, truncates the formulation
    38                                           
    39                                               .. math: Q_{k, p} = (1 - a)^{p / 2} \\sum_{j = 0}^\\infty a^j \\frac{(p / 2)_j}{j!} \\F_{\\chi^2(k + 2 j)}(z + s_{\\min}),
    40                                           
    41                                               where :math:`(x)_j` is the Pochhammer symbol, defined as
    42                                               :math:`(x)_j = x (x + 1) ... (x + j - 1)`, :math:`\\F_k` is the cumulative
    43                                               distribution function of the :math:`\\chi^2(k)` distribution, and
    44                                               :math:`a = s_{\\min} / (z + s_{\\min})`. This is Equation (28) of
    45                                               :cite:p:`hillier2009conditional` or Equation (41) of :cite:p:`hillier2009exact`.
    46                                               The truncation is done such that the error is bounded by a tolerance ``tol``.
    47                                           
    48                                               Uses numerical integration by default.
    49                                           
    50                                               Parameters
    51                                               ----------
    52                                               p: int
    53                                                   Degrees of freedom of the first chi-squared random variable.
    54                                               q: int
    55                                                   Total degrees of freedom.
    56                                               s_min: float
    57                                                   Identification measure.
    58                                               z: float
    59                                                   Test statistic.
    60                                               method: str, optional, default: "numerical_integration"
    61                                                   Method to approximate the critical value function. Must be
    62                                                   ``"numerical_integration"`` or ``"power_series"``.
    63                                               tol: float, optional, default: 1e-6
    64                                                   Tolerance for the approximation of the cdf of the critical value function and
    65                                                   thus the p-value.
    66                                           
    67                                               References
    68                                               ----------
    69                                               .. bibliography::
    70                                                  :filter: False
    71                                           
    72                                                  hillier2009conditional
    73                                                  hillier2009exact
    74                                               """
    75       128         34.0      0.3      0.1      if z <= 0:
    76                                                   return 1
    77                                           
    78       128         20.0      0.2      0.1      if s_min <= 0:
    79                                                   return 1 - scipy.stats.chi2(q).cdf(z)
    80                                           
    81       128         15.0      0.1      0.1      if q < p:
    82                                                   raise ValueError("q must be greater than or equal to p.")
    83       128         14.0      0.1      0.1      if p == q:
    84                                                   return 1 - scipy.stats.chi2(q).cdf(z)
    85                                           
    86       128         22.0      0.2      0.1      if method in ["numerical_integration"]:
    87       128         45.0      0.4      0.2          alpha = (q - p) / 2.0
    88       128         18.0      0.1      0.1          beta = p / 2.0
    89       128         14.0      0.1      0.1          k = q / 2.0
    90       128         29.0      0.2      0.1          z_over_2 = z / 2.0
    91                                           
    92       128         31.0      0.2      0.1          a = s_min / (z + s_min)
    93       128        797.0      6.2      2.9          const = np.power(a, -alpha - beta + 1) / scipy.special.beta(alpha, beta)
    94                                           
    95       128         60.0      0.5      0.2          def integrand(y):
    96                                                       return const * scipy.special.gammainc(k, z_over_2 / y)
    97                                           
    98       256      26033.0    101.7     95.4          res = scipy.integrate.quad(
    99       128         13.0      0.1      0.0              integrand,
   100       128         29.0      0.2      0.1              1 - a,
   101       128         10.0      0.1      0.0              1,
   102       128         14.0      0.1      0.1              weight="alg",
   103       128         29.0      0.2      0.1              wvar=(beta - 1, alpha - 1),
   104       128         14.0      0.1      0.1              epsabs=tol,
   105                                                   )
   106       128         48.0      0.4      0.2          return 1 - res[0]
   107                                           
   108                                               elif method == "power_series":
   109                                                   a = s_min / (z + s_min)
   110                                           
   111                                                   p_value = 0
   112                                           
   113                                                   # Equal to (1 - a)^{p / 2} * a^j * (m/2)_j / j!, where (x)_j is the Pochhammer
   114                                                   # symbol, defined as (x)_j = x (x + 1) ... (x + j - 1). See end of Section 1.0 of
   115                                                   # Hillier's "Exact properties of the conditional likelihood ratio test in an IV
   116                                                   # regression model"
   117                                                   factor = 1
   118                                                   j = 0
   119                                           
   120                                                   # In the Appendix of Hillier's paper, they show that the error when truncating the
   121                                                   # infinite sum at j = J is bounded by a^J * (m/2)_J / J!, which is equal to
   122                                                   # `factor` here. However, their claim
   123                                                   # "F_1(r+1, 1 - m/2, r+2, a) is less than 1 for 0 <= a <= 1"
   124                                                   # is incorrect for m = 1, where F_1(r+1, 1 - m/2, r+2, a) <= 1 / (1 - a) via the
   125                                                   # geometric series. Thus, the error is bounded by a^J * (m/2)_J / J! / (1 - a).
   126                                                   # As G_k(z + l), the c.d.f of a chi^2(k), is decreasing in k, one can
   127                                                   # keep the term G_{k + 2J}(z + l) from the first sum. Thus, we can stop when
   128                                                   # `delta / (1 - a) = factor * G_{k + 2J}(z + l) / (1 - a)` is smaller than the
   129                                                   # desired tolerance.
   130                                                   delta = scipy.stats.chi2(q).cdf(z + s_min)
   131                                                   p_value += delta
   132                                           
   133                                                   sqrt_minus_log_a = np.sqrt(-np.log(a))
   134                                                   tol = tol / (1 + (1 - scipy.special.erf(sqrt_minus_log_a)) / sqrt_minus_log_a)
   135                                           
   136                                                   while delta >= tol:
   137                                                       factor *= (p / 2 + j) / (j + 1) * a
   138                                                       delta = scipy.stats.chi2(q + 2 * j + 2).cdf(z + s_min) * factor
   139                                           
   140                                                       p_value += delta
   141                                           
   142                                                       j += 1
   143                                                       if j > 10000:
   144                                                           raise RuntimeError("Failed to converge.")
   145                                           
   146                                                   p_value *= (1 - a) ** (p / 2)
   147                                           
   148                                                   return 1 - p_value
   149                                           
   150                                               else:
   151                                                   raise ValueError(
   152                                                       "method argument should be 'numerical_integration' or 'power_series'. "
   153                                                       f"Got {method}."
   154                                                   )

Total time: 4.53495 s
File: /Users/mlondschien/code/ivmodels/ivmodels/tests/conditional_likelihood_ratio.py
Function: conditional_likelihood_ratio_test at line 157

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   157                                           @line_profiler.profile
   158                                           def conditional_likelihood_ratio_test(
   159                                               Z,
   160                                               X,
   161                                               y,
   162                                               beta,
   163                                               W=None,
   164                                               C=None,
   165                                               fit_intercept=True,
   166                                               method="numerical_integration",
   167                                               tol=1e-4,
   168                                               critical_values="kleibergen",
   169                                           ):
   170                                               """
   171                                               Perform the conditional likelihood ratio test for ``beta``.
   172                                           
   173                                               If ``W`` is ``None``, the test statistic is defined as
   174                                           
   175                                               .. math::
   176                                           
   177                                                  \\mathrm{CLR}(\\beta) &:= (n - k) \\frac{ \\| P_Z (y - X \\beta) \\|_2^2}{ \\| M_Z (y - X \\beta) \\|_2^2} - (n - k) \\frac{ \\| P_Z (y - X \\hat\\beta_\\mathrm{LIML}) \\|_2^2 }{ \\| M_Z (y - X \\hat\\beta_\\mathrm{LIML}) \\|_2^2 } \\\\
   178                                                  &= k \\ \\mathrm{AR}(\\beta) - k \\ \\min_\\beta \\mathrm{AR}(\\beta),
   179                                           
   180                                               where :math:`P_Z` is the projection matrix onto the column space of :math:`Z`,
   181                                               :math:`M_Z = \\mathrm{Id} - P_Z`, and :math:`\\hat\\beta_\\mathrm{LIML}` is the LIML
   182                                               estimator of :math:`\\beta` (see :py:class:`~ivmodels.KClass`), minimizing the
   183                                               Anderson-Rubin test statistic :math:`\\mathrm{AR}(\\beta)`
   184                                               (see :py:func:`~ivmodels.tests.anderson_rubin_test`) at
   185                                           
   186                                               .. math:: \\mathrm{AR}(\\hat\\beta_\\mathrm{LIML}) = \\frac{n - k}{k} \\lambda_\\mathrm{min}( (X \\ y)^T M_Z (X \\ y))^{-1} (X \\ y)^T P_Z (X \\ y) ).
   187                                           
   188                                               Let
   189                                           
   190                                               .. math:: \\tilde X(\\beta) := X - (y - X \\beta) \\cdot \\frac{(y - X \\beta)^T M_Z X}{(y - X \\beta)^T M_Z (y - X \\beta)}
   191                                           
   192                                               and
   193                                           
   194                                               .. math:: s_\\mathrm{min}(\\beta) := (n - k) \\cdot \\lambda_\\mathrm{min}((\\tilde X(\\beta)^T M_Z \\tilde X(\\beta))^{-1} \\tilde X(\\beta)^T P_Z \\tilde X(\\beta)).
   195                                           
   196                                               Then, conditionally on :math:`s_\\mathrm{min}(\\beta_0)`, the statistic
   197                                               :math:`\\mathrm{CLR(\\beta_0)}` is asymptotically bounded from above by a random
   198                                               variable that is distributed as
   199                                           
   200                                               .. math:: \\frac{1}{2} \\left( Q_{m_X} + Q_{k - m_X} - s_\\mathrm{min} + \\sqrt{ (Q_{m_X} + Q_{k - m_X}  - s_\\mathrm{min})^2 + 4 Q_{m_X} s_\\textrm{min} } \\right),
   201                                           
   202                                               where :math:`Q_{m_X} \\sim \\chi^2(m_X)` and
   203                                               :math:`Q_{k - m_X} \\sim \\chi^2(k - m_X)` are independent chi-squared random
   204                                               variables. This is robust to weak instruments. If identification is strong, that is
   205                                               :math:`s_\\mathrm{min}(\\beta_0) \\to \\infty`, the conditional likelihood ratio
   206                                               test is equivalent to the likelihood ratio test
   207                                               (see :py:func:`~ivmodels.tests.likelihood_ratio_test`), with :math:`\\chi^2(m_X)`
   208                                               limiting distribution.
   209                                               If identification is weak, that is :math:`s_\\mathrm{min}(\\beta_0) \\to 0`, the
   210                                               conditional likelihood ratio test is equivalent to the Anderson-Rubin test
   211                                               (see :py:func:`~ivmodels.tests.anderson_rubin_test`) with :math:`\\chi^2(k)` limiting
   212                                               distribution.
   213                                               See :cite:p:`moreira2003conditional` for details.
   214                                           
   215                                               If ``W`` is not ``None``, the test statistic is defined as
   216                                           
   217                                               .. math::
   218                                                  \\mathrm{CLR(\\beta)} &:= (n - k) \\min_\\gamma \\frac{ \\| P_Z (y - X \\beta - W \\gamma) \\|_2^2}{ \\| M_Z (y - X \\beta - W \\gamma) \\|_2^2} - (n - k) \\min_{\\beta, \\gamma} \\frac{ \\| P_Z (y - X \\beta - W \\gamma) \\|_2^2 }{ \\| M_Z (y - X \\beta - W \\gamma) \\|_2^2 } \\\\
   219                                                  &= (n - k) \\frac{ \\| P_Z (y - X \\beta - W \\hat\\gamma_\\textrm{liml}) \\|_2^2}{ \\| M_Z (y - X \\beta - W \\hat\\gamma_\\textrm{liml}) \\|_2^2} - (n - k) \\frac{ \\| P_Z (y - (X \\ W) \\hat\\delta_\\mathrm{liml}) \\|_2^2 }{ \\| M_Z (y - (X \\ W) \\hat\\delta_\\mathrm{liml}) \\|_2^2 },
   220                                           
   221                                               where :math:`\\hat\\gamma_\\mathrm{LIML}` is the LIML estimator of :math:`\\gamma`
   222                                               (see :py:class:`~ivmodels.KClass`) using instruments :math:`Z`, endogenous
   223                                               covariates :math:`W`, and outcomes :math:`y - X \\beta` and
   224                                               :math:`\\hat\\delta_\\mathrm{LIML}` is the LIML estimator of
   225                                               :math:`(\\beta, \\gamma)` using instruments :math:`Z`, endogenous covariates
   226                                               :math:`(X \\ W)`, and outcomes :math:`y`.
   227                                               Let
   228                                           
   229                                               .. math:: \\Sigma_{X, W, y} := ((X \\ \\ W \\ \\ y)^T M_Z (X \\ \\ W \\ \\ y))^{-1} (X \\ \\ W \\ \\ y)^T P_Z (X \\ \\ W \\ \\ y)
   230                                           
   231                                               and
   232                                           
   233                                               .. math:: \\Sigma_{W, y - X \\beta} := ((W \\ \\ y - X \\beta)^T M_Z (W \\ \\ y - X \\beta))^{-1} (W \\ \\ y - X \\beta)^T P_Z (W \\ \\ y - X \\beta)
   234                                           
   235                                               and
   236                                           
   237                                               .. math:: s_\\mathrm{min}(\\beta) := \\lambda_1(\\Sigma_{X, W, y}) + \\lambda_2(\\Sigma_{X, W, y}) - \\lambda_1(\\Sigma_{W, y - X \\beta}),
   238                                           
   239                                               where :math:`\\lambda_1` and :math:`\\lambda_2` are the smallest and second smallest
   240                                               eigenvalues, respectively.
   241                                               Note that
   242                                               :math:`\\lambda_1(\\Sigma_{X, W, y}) = \\min_{\\beta, \\gamma} \\frac{ \\| P_Z (y - X \\beta - W \\gamma) \\|_2^2 }{ \\| M_Z (y - X \\beta - W \\gamma) \\|_2^2 }`
   243                                               and
   244                                               :math:`\\lambda_1(\\Sigma_{W, y - X \\beta}) = \\min_\\gamma \\frac{ \\| P_Z (y - X \\beta - W \\gamma) \\|_2^2}{ \\| M_Z (y - X \\beta - W \\gamma) \\|_2^2}`.
   245                                           
   246                                               :cite:t:`kleibergen2021efficient` conjectures and motivates that, conditionally on :math:`s_\\mathrm{min}(\\beta_0)`, the statistic
   247                                               :math:`\\mathrm{CLR(\\beta_0)}` is asymptotically bounded from above by a random
   248                                               variable that is distributed as
   249                                           
   250                                               .. math:: \\frac{1}{2} \\left( Q_{m_X} + Q_{k - m_X - m_W} - s_\\mathrm{min}(\\beta_0) + \\sqrt{ (Q_{m_X} + Q_{k - m_X - m_W}  - s_\\mathrm{min}(\\beta_0))^2 + 4 Q_{m_X} s_\\textrm{min} } \\right),
   251                                           
   252                                               where :math:`Q_{m_X} \\sim \\chi^2(m_X)` and
   253                                               :math:`Q_{k - m_X - m_W} \\sim \\chi^2(k - m_X - m_W)` are independent chi-squared random
   254                                               variables. This is robust to weak instruments. If identification is strong, that is
   255                                               :math:`s_\\mathrm{min}(\\beta_0) \\to \\infty`, the conditional likelihood ratio
   256                                               test is equivalent to the likelihood ratio test
   257                                               (see :py:func:`~ivmodels.tests.likelihood_ratio_test`), with :math:`\\chi^2(m_X)`
   258                                               limiting distribution.
   259                                               If identification is weak, that is :math:`s_\\mathrm{min}(\\beta_0) \\to 0`, the
   260                                               conditional likelihood ratio test is equivalent to the Anderson-Rubin test
   261                                               (see :py:func:`~ivmodels.tests.anderson_rubin_test`) with :math:`\\chi^2(k - m_W)`
   262                                               limiting distribution.
   263                                               See :cite:p:`kleibergen2021efficient` for details.
   264                                           
   265                                               Parameters
   266                                               ----------
   267                                               Z: np.ndarray of dimension (n, k)
   268                                                   Instruments.
   269                                               X: np.ndarray of dimension (n, mx)
   270                                                   Regressors.
   271                                               y: np.ndarray of dimension (n,)
   272                                                   Outcomes.
   273                                               beta: np.ndarray of dimension (mx,)
   274                                                   Coefficients to test.
   275                                               W: np.ndarray of dimension (n, mw) or None, optional, default = None
   276                                                   Endogenous regressors not of interest.
   277                                               C: np.ndarray of dimension (n, mc) or None, optional, default = None
   278                                                   Exogenous regressors not of interest.
   279                                               fit_intercept: bool, optional, default: True
   280                                                   Whether to include an intercept. This is equivalent to centering the inputs.
   281                                               method: str, optional, default: "numerical_integration"
   282                                                   Method to approximate the critical value function. Must be
   283                                                   ``"numerical_integration"`` or ``"power_series"``. See
   284                                                   :py:func:`~conditional_likelihood_ratio_critical_value_function`.
   285                                               tol: float, optional, default: 1e-6
   286                                                   Tolerance for the approximation of the cdf of the critical value function and
   287                                                   thus the p-value.
   288                                           
   289                                               Returns
   290                                               -------
   291                                               statistic: float
   292                                                   The test statistic :math:`\\mathrm{CLR}(\\beta)`.
   293                                               p_value: float
   294                                                   The p-value of the test, correct up to tolerance ``tol``.
   295                                           
   296                                               Raises
   297                                               ------
   298                                               ValueError:
   299                                                   If the dimensions of the inputs are incorrect.
   300                                           
   301                                               References
   302                                               ----------
   303                                               .. bibliography::
   304                                                  :filter: False
   305                                           
   306                                                  moreira2003conditional
   307                                                  kleibergen2021efficient
   308                                               """
   309       128        664.0      5.2      0.0      Z, X, y, W, C, beta = _check_test_inputs(Z, X, y, W=W, C=C, beta=beta)
   310                                           
   311       128         31.0      0.2      0.0      n, k = Z.shape
   312       128         35.0      0.3      0.0      mx = X.shape[1]
   313       128         15.0      0.1      0.0      mw = W.shape[1]
   314                                           
   315       128         13.0      0.1      0.0      if fit_intercept:
   316                                                   C = np.hstack([np.ones((n, 1)), C])
   317                                           
   318       128         31.0      0.2      0.0      if C.shape[1] > 0:
   319                                                   X, y, Z, W = oproj(C, X, y, Z, W)
   320                                           
   321                                               # Z = scipy.linalg.qr(Z, mode="economic")[0]
   322       128    2267173.0  17712.3     50.0      X_proj, y_proj, W_proj = proj(Z, X, y, W)  # , orthogonal=True)
   323                                           
   324       128         77.0      0.6      0.0      if mw == 0:
   325                                                   residuals = y - X @ beta
   326                                                   residuals_proj = y_proj - X_proj[:, :mx] @ beta
   327                                                   residuals_orth = residuals - residuals_proj
   328                                           
   329                                                   Sigma = (residuals_orth.T @ X) / (residuals_orth.T @ residuals_orth)
   330                                                   Xt = X - np.outer(residuals, Sigma)
   331                                                   Xt_proj = X_proj - np.outer(residuals_proj, Sigma)
   332                                                   Xt_orth = Xt - Xt_proj
   333                                                   s_min = np.real(
   334                                                       scipy.linalg.eigvalsh(
   335                                                           a=Xt_proj.T @ Xt_proj, b=Xt_orth.T @ Xt_orth, subset_by_index=[0, 0]
   336                                                       )[0]
   337                                                   ) * (n - k - C.shape[1])
   338                                           
   339                                                   # TODO: This can be done with efficient rank-1 updates.
   340                                                   ar_min = KClass.ar_min(X=X, y=y, Z=Z)
   341                                                   ar = residuals_proj.T @ residuals_proj / (residuals_orth.T @ residuals_orth)
   342                                           
   343                                                   statistic = (n - k - C.shape[1]) * (ar - ar_min)
   344                                           
   345       128         25.0      0.2      0.0      elif mw > 0:
   346       128       1175.0      9.2      0.0          XWy = np.concatenate([X, W, y.reshape(-1, 1)], axis=1)
   347       128        507.0      4.0      0.0          XWy_proj = np.concatenate([X_proj, W_proj, y_proj.reshape(-1, 1)], axis=1)
   348                                           
   349       256        605.0      2.4      0.0          XWy_eigenvals = np.sort(
   350       256        265.0      1.0      0.0              np.real(
   351       256       8532.0     33.3      0.2                  scipy.linalg.eigvalsh(
   352       128       1520.0     11.9      0.0                      a=XWy_proj.T @ XWy,
   353       128       1133.0      8.9      0.0                      b=(XWy - XWy_proj).T @ XWy,
   354       128         32.0      0.2      0.0                      subset_by_index=[0, 1],
   355                                                           )
   356                                                       )
   357                                                   )
   358       128    2217473.0  17324.0     48.9          kclass = KClass(kappa="liml").fit(X=W, y=y - X @ beta, Z=Z)
   359       256       5993.0     23.4      0.1          ar = kclass.ar_min(
   360       128        967.0      7.6      0.0              X=W, y=y - X @ beta, X_proj=W_proj, y_proj=y_proj - X_proj @ beta
   361                                                   )
   362                                           
   363       128        125.0      1.0      0.0          statistic = (n - k - C.shape[1]) * (ar - XWy_eigenvals[0])
   364                                           
   365       128         48.0      0.4      0.0          if critical_values == "kleibergen":
   366       128        101.0      0.8      0.0              s_min = (n - k - C.shape[1]) * (XWy_eigenvals[0] + XWy_eigenvals[1] - ar)
   367                                                   else:
   368                                                       XW = np.concatenate([X, W], axis=1)
   369                                                       XW_proj = np.concatenate([X_proj, W_proj], axis=1)
   370                                           
   371                                                       residuals = y - X @ beta - kclass.predict(X=W)
   372                                                       residuals_proj = proj(Z, y_proj - X_proj @ beta - kclass.predict(X=W_proj))
   373                                                       residuals_orth = residuals - residuals_proj
   374                                           
   375                                                       Sigma = (residuals_orth.T @ XW) / (residuals_orth.T @ residuals_orth)
   376                                                       XWt = XW - np.outer(residuals, Sigma)
   377                                                       XWt_proj = XW_proj - np.outer(residuals_proj, Sigma)
   378                                                       s_min = (n - k - C.shape[1]) * min(
   379                                                           np.real(
   380                                                               scipy.linalg.eigvalsh(
   381                                                                   a=XWt_proj.T @ XWt,
   382                                                                   b=(XWt - XWt_proj).T @ XWt,
   383                                                                   subset_by_index=[0, 0],
   384                                                               )
   385                                                           )
   386                                                       )
   387                                           
   388       256      28353.0    110.8      0.6      p_value = conditional_likelihood_ratio_critical_value_function(
   389       128         46.0      0.4      0.0          mx, k - mw, s_min, statistic, method=method, tol=tol
   390                                               )
   391                                           
   392       128         15.0      0.1      0.0      return statistic, p_value

conditional_likelihood_ratio_critical_value_function is less than 1% of runtime.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 93.21%. Comparing base (f5bbf5b) to head (9997f50).
Report is 36 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main      #80      +/-   ##
==========================================
+ Coverage   93.20%   93.21%   +0.01%     
==========================================
  Files          32       32              
  Lines        1604     1607       +3     
==========================================
+ Hits         1495     1498       +3     
  Misses        109      109              

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.