#342 arima, exact ML: --opg (and set arma_vcv op) does not give the outer-product-of-gradient covariance, and --robust inherits it.

[Kirill Naralenkov]

gretl 2026d, Windows 11 Pro 64-bit.

The Command Reference says, under arima, that with native exact ML "the --opg option forces use of the OPG method". The User's Guide (section 26.3) defines that method as (G′G)⁻¹, with G the T × k matrix of per-observation contributions to the gradient, and makes it mle's default. On the same exact log-likelihood, gretl's own mle gives the OPG errors that Stata and statsmodels give, while arima --opg gives other numbers. The point estimates, the log-likelihood and the Hessian-based errors agree to every digit.

The attached arima-opg.inp uses gretl's bundled arma and bjg data:

• it runs arima with its default, --opg, set arma_vcv op and --robust;
• it then runs mle on the per-observation log-likelihood llt that kfilter returns for the same state-space model.

Case A is the User's Guide's Listing 36.1, ARMA(1,1). Case B is the airline model, the example in the arima entry of Stata's Time-Series Reference Manual. The core of case B:

open bjg
arima 0 1 1 ; 0 1 1 ; lg --nc --opg
series w = sdiff(diff(lg))
smpl w --no-missing
# state (e_t, ..., e_t-13): w_t = e_t + th e_t-1 + Th e_t-12 + th Th e_t-13
matrix TT = zeros(14, 14)
loop i = 2..14
    TT[i, i - 1] = 1
endloop
matrix e1 = zeros(14, 1)
e1[1] = 1
bundle ab = ksetup(w, e1, TT, zeros(14, 14))
# air_filter() puts (1, th, 0, ..., 0, Th, th*Th)' into ab.obsymat and sg^2 into ab.statevar[1,1], then calls kfilter
mle logl = ERR ? NA : ab.llt
    ERR = air_filter(&ab, th, Th, sg)
    params th Th sg
end mle

Results, standard errors in parentheses:

A. arma.gdt, arima 1 1 ; y --nc                   phi                   theta
   arima, default (Hessian)                    0.654960 (0.035398)    0.493380 (0.042600)
   mle on kfilter's llt, --hessian             0.654960 (0.035398)    0.493380 (0.042600)
   arima --opg                                 0.654960 (0.034955)    0.493380 (0.040292)
   mle on kfilter's llt, default (OPG)         0.654960 (0.033332)    0.493380 (0.040421)
   lnL -887.9102 in both

B. bjg, arima 0 1 1 ; 0 1 1 ; lg --nc            ma                    sma
   arima, default (Hessian)                   -0.4018231 (0.0896447)  -0.5569365 (0.0731051)
   mle ... --hessian                          -0.4018231 (0.0896447)  -0.5569365 (0.0731051)
   arima --opg, and set arma_vcv op           -0.4018231 (0.0807212)  -0.5569365 (0.0847135)
   mle ... default (OPG)                      -0.4018231 (0.0730299)  -0.5569365 (0.0963116)
   arima --robust                             -0.4018231 (0.0996608)  -0.5569365 (0.0634191)
   mle ... --robust                           -0.4018231 (0.1119815)  -0.5569365 (0.0556247)
   lnL 244.6965 in both, T = 131

Two other programs give mle's OPG, not arima's:

• Stata 19: arima lnair, arima(0,1,1) sarima(0,1,1,12) noconstant, whose default is vce(opg), prints .0730307 and .0963129 (Time-Series Reference Manual, arima, p. 100).
• statsmodels 0.15: SARIMAX(..., cov_type="opg") gives 0.0730298 and 0.0963116 in case B, and 0.033332 and 0.040421 in case A.

--robust is built on the same matrix. arima --robust equals VH·VG⁻¹·VH to every digit (0.0996608, 0.0634191), where VH is the default (Hessian) covariance and VG the --opg one. With the true OPG, mle --robust gives 0.1119815 and 0.0556247.

What --opg computes instead. It is σ̂²(J′J)⁻¹, where J is the T × k matrix of derivatives of the standardized one-step prediction errors v_t/√f_t with respect to the ARMA coefficients (v_t the innovations and f_t their variances from the Kalman filter with σ² = 1, σ̂² the mean of (v_t/√f_t)²). This is a Gauss–Newton matrix, as in least squares, not an outer product of the score contributions. Section C of the attached script computes it with kfilter and fdjac and reproduces arima --opg in both cases: 0.03495546 and 0.04029163, then 0.08072118 and 0.08471348, the largest relative gap in the covariance matrix 3·10⁻⁸ and 2·10⁻⁶. It is the kind of matrix that Box & Jenkins (1976, Time Series Analysis: Forecasting and Control, rev. ed., §7.2.2, p. 233) give as the variance of their least-squares estimates, (X′X)⁻¹σ̂²_a, with the standardized innovations in place of their back-forecast residuals, which give 0.0803 and 0.0774 in case B.

It is asymptotically equivalent to the information matrix. In Harvey (1989, Forecasting, Structural Time Series Models and the Kalman Filter, eq. (3.4.69), p. 142) the information matrix has the terms

½ Σ_t tr[F_t⁻¹ ∂F_t/∂ψ_i F_t⁻¹ ∂F_t/∂ψ_j] + E Σ_t (∂v_t/∂ψ_i)′ F_t⁻¹ (∂v_t/∂ψ_j)

The conditional expectation of J′J/σ² has weight ¼ on the first term, not ½, and statsmodels' cov_type="oim" is the version without the expectation, Harvey's (5.2.2), p. 240. These two differ from each other only through the early observations, where f_t ≠ 1 and the ∂f_t terms are not zero. Like the OPG, both estimate the information matrix, so all three agree in large samples if the model is right, and they differ on a short or strongly seasonal series like bjg.

So the numbers are a legitimate estimate, but not the OPG that the reference names. Either the reference could describe them as a Gauss–Newton (or information-matrix) estimate, or --opg, and --robust with it, should use the score contributions, as mle does. --robust matters more: a sandwich needs the outer product of the scores as its meat, and gretl's uses (J′J)/σ̂².

[Allin Cottrell]

Thanks for your precise and helpful report, Kirill. There's now a fix in gretl git -- taking the option of keeping the --opg label while computing the "genuine" OPG matrix internally rather than what, as you pointed out, is really a Gauss-Newton matrix. After some more testing this change will make its way into the snapshots.

[Allin Cottrell]

Snapshots are now on the way. I'll close this ticket shortly, but if problems with the new code emerge, it can be reopened.