$ontext

  Longley Linear Least Squares benchmark problem.
  Load variance-covariance matrix from gdx file.

  Erwin Kalvelagen, dec 2004


  References:
       http://www.itl.nist.gov/div898/strd/lls/lls.shtml

       Longley, J. W. (1967).
       An Appraisal of Least Squares Programs for the
       Electronic Computer from the Viewpoint of the User.
       Journal of the American Statistical Association, 62, pp. 819-841.

$offtext

set i 'cases' /i1*i16/;
set v 'variables' /empl,const,gnpdefl,gnp,unempl,army,pop,year/;
set indep(v) 'independent variables' /const,gnpdefl,gnp,unempl,army,pop,year/;
set depen(v) 'dependent variables' /empl/;

table data(i,v)
         empl gnpdefl   gnp    unempl    army     pop    year
  i1    60323   83.0   234289   2356     1590    107608  1947
  i2    61122   88.5   259426   2325     1456    108632  1948
  i3    60171   88.2   258054   3682     1616    109773  1949
  i4    61187   89.5   284599   3351     1650    110929  1950
  i5    63221   96.2   328975   2099     3099    112075  1951
  i6    63639   98.1   346999   1932     3594    113270  1952
  i7    64989   99.0   365385   1870     3547    115094  1953
  i8    63761  100.0   363112   3578     3350    116219  1954
  i9    66019  101.2   397469   2904     3048    117388  1955
  i10   67857  104.6   419180   2822     2857    118734  1956
  i11   68169  108.4   442769   2936     2798    120445  1957
  i12   66513  110.8   444546   4681     2637    121950  1958
  i13   68655  112.6   482704   3813     2552    123366  1959
  i14   69564  114.2   502601   3931     2514    125368  1960
  i15   69331  115.7   518173   4806     2572    127852  1961
  i16   70551  116.9   554894   4007     2827    130081  1962
;

data(i,'const') = 1;

alias(indep,j,jj,k);

variables
   b(indep)  'parameters to be estimated'
   sse
;

equation
   fit(i)    'equation to fit'
   sumsq
;

sumsq..   sse =n= 0;
fit(i)..  data(i,'empl') =e= sum(indep, b(indep)*data(i,indep));

option lp = ls;
model leastsq /fit,sumsq/;
solve leastsq using lp minimizing sse;
option decimals=8;
display b.l;


parameter covariance(v,v);

execute_load 'ls.gdx',covariance=covar;
display covariance;