BUG: Fix the layout of the standard errors of LinearFactorModelGMM

linearmodels/asset_pricing/model.py

@@ -1045,17 +1045,21 @@ def fit(
 
         # 2. Step 1 using w = inv(s) from SV
         callback = callback_factory(self._j, args, disp=disp)
-        _default_options: dict[str, Any] = {"callback": callback}
         options = {"maxiter": max_iter}
         opt_options = {} if opt_options is None else opt_options
         options.update(opt_options.get("options", {}))
-        _default_options.update(opt_options)
-        _default_options["options"] = options
+
+        def minimize_kwargs(
+            callback: Callable[[linearmodels.typing.data.Float64Array], None],
+        ) -> dict[str, Any]:
+            # The options of the user apply to every minimization
+            return {"callback": callback, **opt_options, "options": options}
+
         opt_res = minimize(
             fun=self._j,
             x0=np.squeeze(sv),
             args=args,
-            **_default_options,
+            **minimize_kwargs(callback),
         )
         params = opt_res.x
         last_obj = opt_res.fun
@@ -1074,8 +1078,7 @@ def fit(
                     self._j,
                     params,
                     args=args,
-                    callback=callback,
-                    options={"maxiter": max_iter},
+                    **minimize_kwargs(callback),
                 )
                 params = opt_res.x
                 obj = opt_res.fun
@@ -1090,8 +1093,7 @@ def fit(
                 self._j_cue,
                 params,
                 args=cue_args,
-                callback=callback,
-                options={"maxiter": max_iter},
+                **minimize_kwargs(callback),
             )
             params = opt_res.x
 
@@ -1125,7 +1127,11 @@ def fit(
         sel = slice(0, (n * (k + 1)), (k + 1))
         alphas = g.mean(0)[sel, None]
         alpha_vcv = s[sel, sel] / nobs
-        stat = self._j(params, excess_returns, w)
+        if use_cue:
+            # The objective that was minimized updates the weight matrix
+            stat = self._j_cue(params, excess_returns, weight_est_instance)
+        else:
+            stat = self._j(params, excess_returns, w)
         jstat = WaldTestStatistic(
             stat, "All alphas are 0", n - k - nrf, name="J-statistic"
         )

linearmodels/tests/asset_pricing/test_linear_factor_gmm.py

@@ -2,8 +2,10 @@
 from numpy.testing import assert_allclose
 import pandas as pd
 import pytest
+from scipy.optimize import minimize
 
-from linearmodels.asset_pricing.model import LinearFactorModelGMM
+import linearmodels.asset_pricing.model as model_module
+from linearmodels.asset_pricing.model import LinearFactorModel, LinearFactorModelGMM
 from linearmodels.tests.asset_pricing._utility import generate_data, get_all
 
 
@@ -110,6 +112,75 @@ def test_linear_model_gmm_cue_smoke(data):
     get_all(res)
 
 
+def cue_moments(theta, p, f):
+    """Moment conditions of the model with excess returns, written out directly"""
+    nobs, nport = p.shape
+    nfactor = f.shape[1]
+    betas = theta[: nport * nfactor].reshape(nport, nfactor)
+    lam = theta[nport * nfactor : nport * nfactor + nfactor]
+    mu = theta[nport * nfactor + nfactor :]
+    fe = f - mu
+    eps = p - (betas @ lam)[None, :] - fe @ betas.T
+    fc = np.column_stack([np.ones(nobs), f])
+    return np.column_stack([np.hstack([fc * eps[:, [i]] for i in range(nport)]), fe])
+
+
+def cue_objective(theta, p, f):
+    g = cue_moments(theta, p, f)
+    gc = g - g.mean(0)
+    w = np.linalg.inv(gc.T @ gc / g.shape[0])
+    return g.shape[0] * float(g.mean(0) @ w @ g.mean(0))
+
+
+def test_linear_model_gmm_cue_j_statistic():
+    # The J statistic is the minimized value of the continuously updated objective,
+    # which uses the weight matrix at the estimates and not at the starting values
+    rng = np.random.default_rng(3)
+    nobs, nport, nfactor = 150, 8, 2
+    f = 0.05 * (
+        0.2 * rng.standard_normal((nobs, 1)) + rng.standard_normal((nobs, nfactor))
+    )
+    f += np.arange(1, nfactor + 1)[None, :] / 1200
+    betas = rng.chisquare(2, (nfactor, nport)) / 2
+    p = f @ betas + 0.03 * rng.standard_normal((nobs, nport))
+    p += np.arange(nport)[None, :] / nport / 200
+    res = LinearFactorModelGMM(p, f).fit(use_cue=True, disp=0)
+    start = LinearFactorModel(p, f).fit()
+    theta0 = np.r_[
+        np.asarray(start.betas).ravel(), np.asarray(start.risk_premia), f.mean(0)
+    ]
+    best = minimize(
+        cue_objective,
+        theta0,
+        args=(p, f),
+        method="BFGS",
+        options={"gtol": 1e-9, "maxiter": 5000},
+    )
+    assert_allclose(res.j_statistic.stat, best.fun, rtol=1e-4)
+
+
+@pytest.mark.parametrize("kwargs", [{"steps": 3}, {"use_cue": True}])
+def test_linear_model_gmm_opt_options_every_step(monkeypatch, kwargs):
+    # The options of the optimizer apply to every minimization of the estimator
+    data = generate_data(nfactor=1, nportfolio=4, nobs=300, output="numpy")
+    calls = []
+    minimize_original = minimize
+
+    def recording_minimize(*args, **kw):
+        calls.append(kw)
+        return minimize_original(*args, **kw)
+
+    monkeypatch.setattr(model_module, "minimize", recording_minimize)
+    opt_options = {"method": "L-BFGS-B", "tol": 1e-8, "options": {"maxcor": 7}}
+    mod = LinearFactorModelGMM(data.portfolios, data.factors)
+    mod.fit(disp=0, opt_options=opt_options, **kwargs)
+    assert len(calls) >= 2
+    for call in calls:
+        assert call["method"] == "L-BFGS-B"
+        assert call["tol"] == 1e-8
+        assert call["options"] == {"maxiter": 1000, "maxcor": 7}
+
+
 @pytest.mark.parametrize("nfactor", [1, 3])
 @pytest.mark.parametrize("risk_free", [False, True])
 def test_linear_model_gmm_std_errors_layout(nfactor, risk_free):