new;
cls;
library pgraph;
pqgwin many;

/* Step 1: DGP   */ 
b1 = 1;
b2 = 2;
b0 = 0|0;      @ null hypothesis @
k = 2;         @ # of regressors @
sig2 = 4;
Tru_b = b1|b2; @ True parameters @
T = 200;       @ # of observations @
x2 = rndu(T,1)*5;  
emat = rndn(T,1)*sqrt(sig2); @ Error term @
ymat = b1 + x2*b2 + emat;  

/* Step 2: Estimation(OLS) */
X = ones(T,1)~x2;     @ T by 2, Independent variables @
Y = ymat;             @ T by 1, dependent variable @
printr = 1;           @ print out the estimation results if 1 @
{bhat,Yhat,residual,sig2hat,stde,t_val,p_val,F_val,R2,R2_,SC,AIC} = ols1(Y,X,printr);

title("Y and fitted Y");
xy(T,Ymat~Yhat);
title("Error term and Residuals");
xy(T,emat~residual);
end;

proc(12) = ols1(Y,X,printr);
local k,T,bhat,Yhat,ehat,sig2hat,varbhat,stde,
      t_val,R,L,gam,F_val,mY,TSS,RSS,R2,SC,AIC;

  T = rows(Y);
  k = cols(X);

  bhat = inv(X'*X)*X'*Y; @ k by 1, Estimates for b1 and b2  @
  Yhat = X*bhat;         @ T by 1, fitted value @
  ehat = Y - Yhat;       @ T by 1, residuals @
  sig2hat = ehat'ehat/(T-k);     @ Estimates of variance @ 
  varbhat = sig2hat*invpd(X'X);  @ k by k, variance of bhat @
  stde = sqrt(diag(varbhat));    @ k by 1, standard error @
  t_val = (bhat-b0)./stde;       @ k by 1, t values @
  p_val = 2*cdftc(t_val,T-k);      @ k by 1, p value @

  mY = meanc(Y);                 @ mean of the dependent variable @
  TSS  = Y'Y - T*mY^2;
  RSS = ehat'ehat;
  R2 = 1 - RSS/TSS;
  R2_ = 1 - (T-1)*RSS/(TSS*(T-k));
  SC = ln(RSS/T) - k/T*ln(T);
  AIC = ln(RSS/T) - 2*k/T;

/* Test H0: All coefficients are zero */
  R = eye(k);
  L = k;           @ # of Restrictions @
  gam = zeros(k,1);
  F_val = (R*bhat - gam)'invpd(R*invpd(X'X)*R')*(R*bhat-gam)/(L*sig2hat); @ F statistic @


  if printr == 1;
/* Results */
"-------------------------------------------------------------------";
"       Parameter       Estimates          S.E.            t value       p value";
 seqa(1,1,k)~bhat~stde~t_val~p_val; 
"-------------------------------------------------------------------";
"S.E. of regression " sqrt(sig2hat);
"F value            " F_val;
"p value           " cdffc(F_val,L,T-k);
"R2            " R2;
"adjusted R2   " R2_;
"SC           " SC;
"AIC          " AIC;
"-------------------------------------------------------------------";
  endif;

  retp(bhat,Yhat,ehat,sig2hat,stde,t_val,p_val,F_val,R2,R2_,SC,AIC);
endp;