new;
cls;

/* Data loading */
t=144; @ 1960:1-2001:4 @
load data[t,3]=c:\data\watson.txt; 

FFR = data[.,1]; @ Federal Funds rate @
inf = data[.,2]; @ inflation @
un = data[.,3];  @ Unemployment rate @

/* Estimatin */
/* Step 1: Unrestricted OLS */
bp = 80;
D = zeros(80,1)|ones(t-80,1);
Y = FFR;
X1 = ones(t,1);
X2 = X1.*D;
X3 = inf;
X4 = inf.*D;
X5 = un;
X6 = un.*D;
X = X1~X2~X3~X4~X5~X6;

printr = 0;

{bhat,Yhat,resid,sig2hat,stde,t_val,p_val,F_val,R2,R2_,SC,AIC} = ols1(Y,X,printr);

/* method 1 : Wald Test */
@  1st method @
R = 0~1~0~0~0~0|
  0~0~0~1~0~0|
  0~0~0~0~0~1;
L = rows(R); @ # of Restrictions @
gam=0|0|0;
F1_val=(R*bhat - gam)'inv(R*inv(X'X)*R')*(R*bhat-gam)/(L*sig2hat); @ F statistic @


/* Step 2: Restricted OLS */
Y = FFR;
X0 = X1~X3~X5;
printr = 0;
{bhat0,Yhat0,resid0,sig2hat0,stde0,t_val0,p_val0,F_val0,R20,R20_,SC0,AIC0} = ols1(Y,X0,printr);


/* method 2: Wald Test */
RSS_res=resid0'resid0; @ Restricted RSS @
RSS_unr=resid'resid; @ Unestricted RSS @
k = cols(X); @ # of regressors in the unrestricted reg @
F2_val=(T-k)*(RSS_res - RSS_unr) / (L*RSS_unr);  @ F statistic @

"-----------------------";
" Result for Chow Test  ";
"-----------------------";
"F value from method 1  " F1_val;
"F value from method 2  " F2_val;
"p value " cdffc(F1_val,L,T-k);

end;



proc(12) = ols1(Y,X,printr);
local k,T,bhat,Yhat,ehat,sig2hat,varbhat,stde,p_val,R2_,
      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./stde;       @ k by 1, t values @
  p_val = 2*cdftc(abs(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;