set verbose off

function matrix logrel_err (const matrix est, const matrix true)
  return -log10(abs(est - true) ./ abs(true))
end function

open longley.gdt -q
list LX = prdefl .. year
ols employ 0 LX -q
matrix b_gretl = $coeff

mwrite({employ} ~ {const} ~ {LX}, "alldata.mat", 1)

foreign language=python
   import numpy as np
   X = gretl_loadmat('alldata.mat', 1)
   # NumPy's OLS   
   b = np.linalg.lstsq(X[:,1:], X[:,0])[0]
   gretl_export(np.transpose(np.matrix(b)), 'py_b.mat', 1)
end foreign

# NIST's certified coefficient values
matrix b_nist = {-3482258.63459582, 15.0618722713733,
    -0.358191792925910E-01, -2.02022980381683,
    -1.03322686717359, -0.511041056535807E-01,
     1829.15146461355}'

matrix b_numpy = mread("py_b.mat", 1)
matrix E = logrel_err(b_gretl, b_nist) ~ logrel_err(b_numpy, b_nist)
cnameset(E, "gretl python")
printf "Log-relative errors, Longley coefficients:\n\n%#12.5g\n", E
printf "Column means\n%#12.5g\n", meanc(E)