$ontext

   Ordinary Least Squares (OLS) by solving
   the normal equations.

   Erwin Kalvelagen, october 2000

$offtext


set i /i1*i40/;

$include expdata.inc

set j 'parameters to be estimated' /
      constant  'constant term'
      coeff1    'income coefficient'
/;
alias (j,jj);


parameters
   X(i,j)   'the X matrix in standard OLS notation (dependent variables)'
   y(i)     'the y vector in standard OLS notation (independent variables)'
   XX(j,jj) "the matrix (X'X)"
   Xy(j)    "the vector (X'y)"
;

X(i,'constant') = 1;
X(i,'coeff1')   = data(i,'income');

y(i) = data(i,'expenditure');

XX(j,jj) = sum(i, X(i,j)*X(i,jj));

Xy(j) = sum(i, X(i,j)*y(i));

equations
    dummy_eq   'dummy objective equation'
    normal(j)  "normal equations (X'X)b = X'y"
;

variables
    b(j)       'parameters to estimate'
    dummy_var  'dummy objective variable'
;

dummy_eq..  dummy_var =e= b('constant');
normal(j).. sum(jj, XX(j,jj)*b(jj)) =e= Xy(j);

model ols2 /dummy_eq,normal/;
solve ols2 using lp minimizing dummy_var;

display b.l;