load datall[421,6]=caschool1.txt;
dat = datall[2:101,.];
vnames = datall[1,.]; 
$vnames;
clear datall;

y = dat[.,1];
t = rows(dat);

ind1 = { 4, 5, 6 } ;
ind2 = { 2 3 } ;
x1 = ones(t,1)~dat[.,ind1];
x2 = dat[.,ind2];
x = x1~x2; 

{b,ehat}=ls(y,x);
{b1,ehat1} = ls(y,x1);

px = x*invpd(x'*x)*x';
p1 = x1*invpd(x1'*x1)*x1';
p2 = x2*invpd(x2'*x2)*x2';
ht = diag(p1);
alpha = ehat1./(1-ht);
" leverage    and   Influence "  ht~alpha;
del = invpd(x1'*x1)*x1';
del[.,1]*alpha[1];

m1 = eye(t) - p1;
m2 = eye(t) - p2;

ys1 = m2*y;
xs1 = m2*x1;

{bs1,ehats1} = ls(ys1,xs1);
 
proc influence(y,x,i);
local mx,px,ehat,ht,alpha,del,d; 
  px = x*invpd(x'*x)*x';
  mx = eye(rows(y))-px;
  ehat = mx*y;  
  ht = diag(px);
  alpha = ehat./(1-ht);
  del = invpd(x1'*x1)*x1';
  d=del[.,i]*alpha[i];
retp(d);

endp;
/* R-square */
proc rsquare(b,x,y);
local t, yhat, ybar, ssr, sst, sse, r2;
  t=rows(y);
  yhat=x*b;
  ybar=meanc(y).*ones(t,1);
  SSR=(yhat-ybar)'*(yhat-ybar);
  SST=(y-ybar)'*(y-ybar);
  SSE=(y-x*b)'*(y-x*b);
  r2=ssr/sst;
  "SSR         SST            SSE            r2 " ;
  ssr~sst~sse~R2;
  retp(r2);
endp;


/*  Wald Test */
proc Wald(y,x,R,lr);
local t,k,j,t_k,b,sig2,w,pval;
  T=rows(X);
  K=cols(X);
  J=rows(R);
  t_k=t-k;
  b=invpd(x'*x)*x'*y;
  sig2=(y-x*b)'*(y-x*b)/(T_K);
  W=(R*b-lr)'*invpd(R*invpd(x'*x)*R')*(R*b-lr)/(J*sig2);
  pval=1-cdffc(w,j,t_k);
  "W     "   "J     "    "DF          "   "pval     "; 
  w~J~T_K~pval;
  retp(w);
endp;

proc RESET(y,x);
local k,b,yhat,e2,e3,z1,r,lr,w1,w2;
  k=cols(x);
  b=invpd(x'*x)*x'*y;
  yhat=x*b;
  e2=yhat.*yhat;
  e3=yhat.*yhat.*yhat;
  z1=x~e2;
    r=zeros(1,k)~1;
    lr=0;
    print "Reset Power 2";
    w1=wald(y,z1,r,lr);
    z1=x~e2~e3;
    r=zeros(2,k)~eye(2);
    lr=0;
    "Reset Power 3";
    w2=wald(y,z1,r,lr);
retp(w1~w2);
endp;


 proc (2) = ls(y,x);
 local t,k,b,sig2,cov,se,tvals,pvals,mx,px; 
  T=rows(X);
  K=cols(X);
  px = x*invpd(x'*x)*x';
  mx = eye(t) - px;
  b=invpd(x'*x)*x'*y;

  /*  Sigma squared */
  sig2=(y-x*b)'*(y-x*b)/(T-K);
  print sig2;

  /* Estimated Covariance Matrix */
  cov=sig2*invpd(x'*x);
  
  /* Standard Errors */
  se=sqrt(diag(cov));
  /* Tvalues */
  tvals=b./se;
  pvals=2*(cdftc(abs(tvals),t-k));

  print "Ordinary Least Squares";
  b~se~tvals~pvals;
retp(b,mx*y);  
endp;

/* a */
{ba,ehata}=ls(y,x2);
sumc(ehat-ehata);

/* b */
{bb,ehatb}=ls(p1*y,x2);
sumc(ehat-ehatb);

/* c */
{bc,ehatc}=ls(p1*y,p1*x2);
sumc(ehat-ehatc);

/* d */
{bd,ehatd}=ls(px*y,x);
sumc(ehat-ehatd);

/* e */
{be,ehate}=ls(px*y,x2);
sumc(ehat-ehate);

/* f */
{bf,ehatf}=ls(m1*y,x2);
sumc(ehat-ehatf);

/* g */
{bg,ehatg}=ls(m1*y,m1*x2);
sumc(ehat-ehatg);

/* h */
{bh,ehath}=ls(m1*y,x1~m1*x2);
sumc(ehat-ehath);

/* i */
{bi,ehati}=ls(m1*y,m1*x2);
sumc(ehat-ehati);

/* j */
{bj,ehatj}=ls(px*y,m1*x2);
sumc(ehat-ehatj);

/* exam a */
{ba,ehata}=ls(px*y,px*x1~px*x2);
sumc(ehat~ehata);