new;
cls;
library pgraph;
pqgwin many;
/* Step 1: DGP   */ 
b1 = 1;
b2 = 2;
b0 = 1|2;  @ 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 @

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 @

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;

/* Results */
"-------------------------------------------------------------------";
"       True values      Estimates            S.E.        t value   ";
Tru_b~bhat~stde~t_val; 
"-------------------------------------------------------------------";
"S.E. of regression is  " sqrt(sig2hat);
"R2  is          " R2;
"adjusted R2  is " R2_;
"SC  is          " SC;
"AIC is          " AIC;
"-------------------------------------------------------------------";

title("Y and fitted Y");
xy(T,Ymat~Yhat);

title("Error term and Residuals");
xy(T,emat~ehat);

end;