%let location=C:\Documents and Settings\woodph\My Documents\glmwood\kutnersolutions\Chapter5;
*For doing regression in matrix form, it's instructive to use proc IML. ;
proc iml;
x={16000 23,
   33000 47,
   21000 35};
print "Refer to page 176" x;
a={1 0, 5 10};print "Refer to page 177" a;
b={4 7 12 16,
   3 15 9 8};print "Refer to page 177" b;

*data sets can also be read directly into IML when it's convenient to use 
 matrix algebra to do the work;
data adata;
input x1 x2;
cards;
2 5 
7 10
3 4
;
proc iml;
use adata;
read all into a;
print "data read in " a;
at=a`;
print "A transpose is" at;
title "transpose of matrix. Refer to page 178";run;
/*You can use proc transpose to transpose a data set if that's all you 
  want to do:*/
proc transpose data=adata out=tadata;
proc print data=adata;
title "original data";
proc print data=tadata;
title "transposed data. Refer to page 178";run;
proc iml;
a={1 4,
   2 5,
   3 6};
b={1 2,
   2 3,
   3 4};
sum=a+b;
diff=a-b;
print "For " a "and" b "sum is" sum;
print "For " a "and" b "difference is" diff;
title "Matrix sum. Refer to page 180";run;
proc iml;
a={2 7, 9 3};
sa=4#a;*scalar multilier is # not *, which is reserved for matrix multiplication;
print "For" a "scalar multiplied by 4 result is" sa;
b={9 25,15 18};dsb=b/3;
print b "divided by scalar 3 is" dsb;
title "multipication by a scalar. Refer to page 182";
proc iml;
a={2 5,4 1};b={4 6,5 8};
ab=a*b;
print a "multiplied by " b "yields" ab;
ba=b*a;
print "ab is not ba necessarily and may not even be defined"
 ab ba;
title "refer to page 184";
proc iml;
X={2, 3, 5};
ss=x`*x;
print "refer to page 184" x ss;
proc iml;
a={1 4 6,
   4 2 5,
   6 5 3};
ta=a`;
print "a is" a "and transpose of a is" ta;
title "symmetric matrix. Refer to page 185";
proc iml;
bv={4,
    1,
	10,
	5};
diagb=diag(bv);
print "vector version of b is:" bv "and diagonal version of bv is" diagb;
square={1 2 3,
        4 5 6,
		7 8 9};
dvec=diag(square);
print "and, if given a square input matrix, diag strips out the diagonal elements e.g.:"
"diag of " square "yields" dvec;
title "diagonal of a mtrix. refer to page 185";run;
proc iml;
x={1 2 3,
   4 5 6,
   7 8 9};
ident=i(3);
xi=x*ident;
ix=ident*x;
print "if x is" x "then ix is" ix "and xi is" xi;
scalar=5#ident;
print "making your own scalar matrix if scalar is 5" scalar;
title "identity matrix. Refer to page 186";
proc iml;
colvec=j(3,1,1);
print "column vector" colvec;
jmat=j(3,3,1);
print "j matrix" jmat;
n=colvec`*colvec;
print "sample size" n;
jmult=colvec*colvec`;
print "j via multiplication" jmult;
zero=j(3,1,0);
print "zero vector" zero;
title "column vectors, j, and sample size calculations. Refer to page 187";
proc iml;
a={1 2 5 1,2 2 10 6,3 4 15 1};
ta=a`;
ca=echelon(ta);*iml tests for row, as opposed to column echelon;
print a ra;
title "rank of a matrix. refer to page 188";
proc iml;
a={2 4,3 1};
ia=inv(a);
print "for " a "inverse is" ia;
v={3 4 2};
da=diag(v);
ida=inv(da);
print "for a" da "inverse is" ida;
title "inverse of a matrix. Refer to page 189-191";
proc iml;
x={2 4,3 1};
const={20,10};
y=inv(x)*const;
print "for " x "*y=" const "solution is" y;
title "solution of a system of linear equations. Refer to page 192";run;
/*p. 196 Bastardizing the code a little bit, we have a discussion of the computation of the variance of a sum
  and difference of two variables and the covariance between the two. To keep the notation consistent
  with the rest of the book, it's easier to consider the matrices in their transposed form.For example:*/
proc iml;
y={-4 -4,
   -3 -3,
   -2 -2,
   -1 -1,
   0 0,
   1 1,
   2 2,
   3 4,
   4 3};
sigmay=y`*y/9;
a={1 1,
   -1 1};
w=y*a;
print "differences and sums" w "of" y;
print "raw variance/covariance matrix" sigmay;
sigmaad=a*sigmay*a`;
print "variance/covariance matrix of differences and sums" sigmaad;
title "Variances of sum and difference scores. Refer to page 196";
proc iml;
PROC IMPORT OUT= WORK.Toluca 
            DATAFILE= "&location\chapter5.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH01TA01$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc iml;
use toluca;
read all into data[colname=vars];
print "raw data is" data "and variable names are" vars;
ones=j(nrow(data),1,1);
namepvar=vars[1,1]||"One";
x=data[,1]||ones;
print "adjoined score matrix" x[colname=namepvar];
xpx=x`*x;

print "P. 200: x`x is" xpx[rowname=namepvar colname=namepvar];
y=data[,2];
xpy=x`*y;
namey=vars[1,2];
print "x`y is" xpy[rowname=namepvar colname=namey];
inv_xpx=inv(x`*x);
print "Page 201: Inverse of x`x is" inv_xpx[rowname=namepvar colname=namepvar];
beta=inv_xpx*xpy;
print "unstandardized regression weight and intercept is" beta[rowname=namepvar colname=namepvar];
beta=inv(x`*x)*x`*y;
print "also calculated in one statement as:" beta[rowname=namepvar colname=namepvar];
ypred=x*beta;
print "Page 202: predicted values of Y" ypred[colname=namey];
hat=x*inv(x`*x)*x`;
print "hat matrix is" hat;
ypred=hat*y;*"yhat is equal to hat y"-a wee bit of humor from your statistician friends :-) ;
print "y pred using hat matrix" ypred;
hh=hat*hat;*no mattter how many times you wear it, it's still a hat;
print "Page 203: For hat of" hat "hh=" hh "(i.e., is ""idempotent"")";
resid=y-x*beta;
print "residual scores are:" resid;
ni=i(nrow(y));
resid=(ni-hat)*y;
print "residual scores via hat matrix:" resid;
mse=(resid`*resid)/(nrow(y)-2);
print "mse is" mse;
vare=mse*(ni*hat);
print "variance/covariance matrix of e's:" vare;
ssy=y`*y;
print "page 204-205: sum of squares of (raw) y" ssy;
j=j(nrow(y),nrow(y),1);
sstot=y`*y-(y`*j*y)/nrow(y);
print "sum of squares total" sstot;
sserror=(y-x*beta)`*(y-x*beta);
print "sum of squares error" sserror;
sserror=y`*y-beta`*x`*y;
print "sum of squares error using 5.84a" sserror;
*ss calculated based on quadratic form (p. 206);
sstotal=y`*(ni-j/nrow(y))*y;
sserror=y`*(ni-hat)*y;
ssreg=y`*(hat-j/nrow(y))*y;
print "quadratic forms of SS yield:" sstotal sserror ssreg;
varb=mse#inv(x`*x);
print "Page 207: Variance/covariance matrix of b's" varb;
*mean response, page 208;
xh={65 1};
varmh=xh*varb*xh`;
print "page 208, variance of mean" varmh;
varpred=mse*(1+xh*inv(x`*x)*xh`);
print "page 209, variance of new observation" varpred;
title "regression of adjoined score matrix. refer to pp. 200-209";