$ontext

   OLS plus statistics

   Erwin Kalvelagen, november 2000

$offtext

set i /i1*i40/;

$include expdata.inc

set j /'constant','coeff1'/;
alias (j,jj);
alias (i,ii);

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= 0;
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 "----------------estimates--------------------",
        b.l;

parameters
   residual(i)   'residuals (errors)'
   yhat(i)       'predicted y'
   A(i,i)        "Theil's A matrix: I - (1/n) (iota*iota')"
;


scalars
   sse           'sum of squared errors'
   sst           'total sum of squares'
   ssr           'regression sum of squares'
   n             'number of observations'
   df            'degrees of freedom'
   r2            'r-square'
   r2adj         'r-square adjusted'
   llf           'log of the likelihood function'
   sigma_squared 'variance of the estimate'
   sigma         'standard error of the estimate'
   pi            '3.1415...'
;

yhat(i) = sum(j, x(i,j)*b.l(j));
residual(i) = y(i) - yhat(i);

sse = sum(i, sqr(residual(i)));

n = card(i);
df = n - 2;
sigma_squared = sse/df;
sigma = sqrt(sigma_squared);
pi = 4*arctan(1);

A(i,ii) = 1$sameas(i,ii) - (1/n);

sst = sum((i,ii), data(i,'expenditure')*A(i,ii)*data(ii,'expenditure'));
ssr = sst - sse;

r2 = 1-sse/sst;

r2adj = 1 - ((n-1)/df)*(1-r2);

llf = -(n/2)*log(2*pi*sse/n)-(n/2);



display "--------------statistics---------------------",
        sse,
        sigma_squared,
        sigma,
        r2,
        r2adj,
        llf;


alias (j,jjj);
variable var(j,jj) 'variance of the estimators';
equation variance(j,jj);
variance(j,jj).. sum(jjj, XX(j,jjj)*var(jjj,jj)) =e= sigma_squared$sameas(j,jj);

model mvar  /dummy_eq,variance/;
solve mvar using lp minimizing dummy_var;

parameters
    se(j)      'standard error'
    t(j)       "t ratio's"
    partial(j)
;

se(j) = sqrt(var.l(j,j));

t(j) =  b.l(j)/se(j);

partial(j) = t(j)/sqrt(sqr(t(j))+df);

display "--------------------------------------------------",
        se,
        t,
        partial
;