This is the program file for Kevin L. Kliesen and Frank A. Schmid's 
   article "Inflation Expectations and Macro Data Releases"
   in the Federal Reserve Bank of St. Louis's Review, May/June 2004 issue.
This file contains 2 Gauss programs: neweywest.prg and amemiya.prg



(neweywest.prg)
@ Gauss program to correct test statistics and perform significance tests.@
@ For a description of the program's functions see the program and comments below @


new; #lineson;
library pgraph;
output file=c:\prg\strips\neweywest.out reset;
str=datestr(0);
tr=timestr(0);

"GAUSS                            time  ";; tr;; "  date (US) ";; str;

?;"input file: neweywest.prg";

_iv=0; @ 1: instrumental variables @
rob=0; @ 1: delete days on which a speech was given (not used in paper); 2: use core measures only @
match=0; @ shorten window for STRIPS/TIIS measure to make it comply with window of CM measure @

ig=0;
do while ig.<2;
ig=ig+1;

if match.==1 and ig.==2; break; endif; @ do STRIPS/TIIS measure only @

white=1;
@ see below for manual activation of newest @

vnames="xx";
dtset0="c:\\data\\strips\\gauss\\newey\\newey";
dtset0="c:\\data\\yield\\gauss\\regress";
declare bvar2;

open f0=^dtset0;
dat=readr(f0,rowsf(f0));
f0=close(f0);
clear f0;

@ obs are currently stacked like this: 1997|1998|...|2003 @

@ ----- define variables ----- @

if ig.==1;
  ?;?;"STRIPS (OTR) INFLATION MEASURE";?;
  if match./=1;
    dat=delif(dat,(dat[.,1].<=1).*(dat[.,2].<31).*(dat[.,3].==1997));
  elseif match.==1; @ use criterion used when ig.==2 @
    dat=delif(dat,dat[.,3].<1999);
    dat=delif(dat,(dat[.,1].<=1).*(dat[.,2].<5).*(dat[.,3].==1999));
  endif;
  y__=dat[.,4]; @ STRIPS inflation measure @
  /* format /rd 10,5; y__[1:10]; */
  dtset1="c:\\data\\strips\\gauss\\kernel\\res_otr";
/*
  dat[.,44]=missrv(dat[.,44],-9999);
  ex=dat[.,44]./=-9999;
  sumc(ex); wait;
*/
elseif ig.==2;
  ?;?;"CONSTANT MATURITY INFLATION MEASURE";?;?;
  dat=delif(dat,dat[.,3].<1999);
  dat=delif(dat,(dat[.,1].<=1).*(dat[.,2].<5).*(dat[.,3].==1999));
  y__=dat[.,5]; @ constant-maturity inflation measure @
  /* format /rd 10,5; y__[1:10]; */
  dtset1="c:\\data\\strips\\gauss\\kernel\\res_cm";
/*
  dat[.,44]=missrv(dat[.,44],-9999);
  ex=dat[.,44]./=-9999;
  sumc(ex); wait;
*/
endif;

@ compile dataset @

@ dat=dat[.,1:43 45 46]; @

@ use all variables bar the trading day indicator @
x__=dat[.,7:cols(dat)]~ones(rows(dat),1); @ standardized macro announcement variables~const @
let seqa_=1 2 3 6 5 8 9 7 10 11 12 13 33 35 34 15 17 16 18 25 19 20 21 23 24 26 29
          30 31 32 27 28 36 37 38 39 4 22 14 40 42;

@ cols(x__); rows(seqa_); wait; @

@ do not use the core numbers @
  x__=dat[.,7:12 14:35 37:39 41:cols(dat)]~ones(rows(dat),1);
                                     @ 12/13: CPI/Core; 35/36: PPI/Core;  39/40: Sales/ex.Auto @
  let seqa_=1 2 3 6 5 8   7 10 11 12 13 33 35 34 15 17 16 18 25 19 20 21 23 24 26 29
          30 31    27 28 36    38 39 4 22 14 40 42;
                                     @ 8/9: CPI/Core; 31/32: PPI/Core;  36/37: Sales/ex.Auto @

if rob.==2; @ use the core numbers @
  x__=dat[.,7:11 13:34 36:38 40:cols(dat)]~ones(rows(dat),1);
                                     @ 12/13: CPI/Core; 35/36: PPI/Core;  39/40: Sales/ex.Auto @
  let seqa_=1 2 3 6 5   9 7 10 11 12 13 33 35 34 15 17 16 18 25 19 20 21 23 24 26 29
          30    32 27 28    37 38 39 4 22 14 40 42;
                                     @ 8/9: CPI/Core; 31/32: PPI/Core;  36/37: Sales/ex.Auto @
endif;

x__=missrv(x__,0); @ turn empty cells into zeros @

y_=rev(y__);
x_=rev(x__);

@ obs are now stacked like this: 2003|2002|...|1997 @

@ now isolate missing observations (there was trading but we do not have the observation) @
@ note that days on which there is no trading have been eliminated from the dataset @

seqas=seqa(1,1,rows(dat)); @ index that includes missing observations @

y_=missrv(y_,-9);
ymiss=indexcat(y_,-9); @ positions of missing values @

y=delif(y_,y_.==-9);
x=delif(x_,y_.==-9);
seqam=delif(seqas,y_.==-9);

if rob.==1;
  y    =delif(y,x[.,cols(x)-1].==1);
  seqam=delif(seqam,x[.,cols(x)-1].==1);
  ymiss=ymiss|indexcat(x[.,cols(x)-1],1);
  x    =delif(x,x[.,cols(x)-1].==1);
endif;

wks=rows(seqas); @ number of obs in time, including missing values @
@ wks=rows(y); @
l=floor(4.*((wks./100).^(2./9))); @ lag for Ljung-Box and Newey-West @

/* ************** deal with speech/testimony indicator ************* */
/*
z=        x[.,cols(x)-2:cols(x)]; @ do not exlude speech/testimony indicator: next to last column @
seqaz=seqa_[cols(x)-2:cols(x)];
zex=(z[.,1].==0).*(z[.,1].==0);
@ z=delif(z,zex); @
@ create intercept indicator variable @
z=z[.,1 2]~zex~z[.,3]; @ two regressors, 0/1 intercept indicator, intercept @
seqaz=seqaz[1 2]|41|seqaz[3];
*/
x_=     x[.,1:cols(x)-2 cols(x)]; @ exlude speech/testimony indicator: next to last column @
seqax=seqa_[1:cols(x)-2 cols(x)];
@ create intercept indicator variable @
xex=(stdc(x_[.,1:cols(x_)-1]').==0);
@ xex=(xex.==0); @
x_=x_[.,1:cols(x_)-1]~xex~x_[.,cols(x_)];
seqax=seqax[1:rows(seqax)-1]|41|seqax[rows(seqax)];
@ "number of nonzero observations:" rows(xex)-sumc(xex); wait; @

/* ************************** core program ************************* */

x=x_;

@ frequency distribution of conincidence in surprises @
/*
e=x[.,1:cols(x)-2];
e_=(e./=0);
sumc(e_');
e__=(sumc(e_'));
format /rd 8,0;
sumc((e__.==0));
sumc((e__.==1));
sumc((e__.==2));
sumc((e__.==3));
sumc((e__.==4));
sumc((e__.==5));
sumc((e__.==6));
sumc((e__.==7));
sumc((e__.==8));
sumc((e__.==9));
wait;
*/
@ for iv approach @
z=x_;
posmed=median(selif(z[.,cols(x_)-2],z[.,cols(x_)-2].>0));
negmed=median(selif(z[.,cols(x_)-2],z[.,cols(x_)-2].<0));
@ posmed; negmed; wait; @

zpos_=z;
zneg_=z;
zpos_[.,cols(x_)-2]= (z[.,cols(x_)-2].>posmed);
zneg_[.,cols(x_)-2]=-(z[.,cols(x_)-2].<negmed);
z[.,cols(x_)-2]=zpos_[.,cols(x_)-2]+zneg_[.,cols(x_)-2];

@ z[.,cols(x_)-2]; wait; @

b=inv(x'x)*x'y;

if _iv.==1;
  b=inv(z'x)*z'*y;
endif;

res=y-x*b; @ residuals @
sige=(res'res)./(rows(x)-cols(x)); @ Varianz der Schaetzung @
bvar=inv(x'x).*sige; @ Varianz der Schaetzer b @

if _iv.==1;
  bvar=sige.*inv(z'x)*(z'*z)*inv(x'z);
endif;

@ save residuals for kernel estimation @

create f1=^dtset1 with ^vnames,1,8;
call writer(f1,res);
f1=close(f1);

@--------------------------------@

@ expand the residual and the x vector @

resnew=res~seqam;
xnew=x~seqam;
resadd=zeros(rows(ymiss),1)~ymiss;
xadd  =zeros(rows(ymiss),cols(x))~ymiss;

res_=resnew|resadd;
x_  =xnew|xadd;

if unique(res_[.,2],1) ne rows(res_); "error (x)--wait"; wait; endif;

res_=sortc(res_,cols(res_));
res_=res_[.,1];

x_=sortc(x_,cols(x_));
x_=x_[.,1:cols(x_)-1];

@ Ljung-Box @

rj=zeros(l,2); @ fuer Ljung-Box; Greene, 2nd, pp. 426-7 @
/*
res_=res;
x_=x;
*/
i2=0; @ Greene: j @
do while i2<l;
i2=i2+1;

rj[i2,2]=wks-i2;
w=1-i2./(l+1);

i3=i2; @ Greene: t @
do while i3<wks;
i3=i3+1;

@ fuer Ljung-Box; Greene, pages 426-27 @
rj[i2,1]=rj[i2,1]+(res_[i3].*res_[i3-i2]);

mat   =  w.*res_[i3].*res_[i3-i2]
       .*(x_[i3,.]'*x_[i3-i2,.]+x_[i3-i2,.]'*x_[i3,.]);
/* i2;; i3; wait; */
if i2.==1 and i3.==i2+1; s_star=mat; else; s_star=s_star+mat; endif;

endo; @ i3 @
endo; @ i2 @

@ calculate Ljung-Box @

rj[.,1]=rj[.,1]./sumc(res.^2); @ rj, wie bei Greene, page 427, erste Formel @
rj=(rj[.,1].^2)./rj[.,2]; @ quadriert und durch T-j dividiert @
qprime=sumc(rj).*(wks).*(wks+2);
p1lj=cdfchic(qprime,l);

if p1lj<0.05; newest=1; else; newest=0; endif; @ serielle Korrelation @

format /rd 5,0; "Ljung-Box-statistic: (";;l;;")";;
format /re 12,3; qprime;; format /rd 8,4; ",  rhs-probability:" p1lj;;
if p1lj le 0.01; "***"; elseif p1lj gt 0.01 and p1lj le 0.05; "**";
elseif p1lj gt 0.05 and p1lj le 0.1; "*"; else; " "; endif;

newest=1;

@--------------------------------@

@ White @

if _iv.==1; x_=x; x=z; endif;

if white.==1;
hat=diagrv(zeros(rows(x),rows(x)),res.^2); @ diagonale Matrix @
/*
bvar2=inv(x'x)*x'*hat*x*inv(x'x); @ White @
*/
@ Vereinfache Matrix-Multiplikation @
@ Greene, 2nd ed., page 391 @
xt=x';
s0 =zeros(cols(x),cols(x));
s0_=zeros(cols(x),rows(x)); @ Zwischenmatrix @
i=0;
do until i==cols(x); @ ueber alle Regressoren @
i=i+1;
j=0;
do until j==rows(x); @ ueber alle Beobachtungen @
j=j+1;
s0_[i,j]=(xt[i,j].*hat[j,j]);
endo; @ ueber j; d.h. Zeile i ist nun vollstaendig @
s0[i,.]=s0_[i,.]*x;
endo; @ ueber i @
endif; @ zu White @

@--------------------------------@

if white .==1; s_ast=s0;        endif;
if newest.==1; s_ast=s0+s_star; endif; @ Greene, page 422 @
@ clear s0,s_star; @

if _iv.==1; x=x_; endif;

@ Greene, 2nd ed., page 408 @
if white.==1 or newest.==1;
bvar2=inv(x'x)*s_ast*inv(x'x);
  if _iv.==1;
  bvar2=inv(z'x)*s_ast*inv(x'z);
  endif;
endif;

@ R-squared @
r2=(b'x'y-rows(y).*(meanc(y).^2))./(y'y-rows(y).*(meanc(y).^2));
r2adj=1-(1-r2) .* (rows(y)-1) ./ (rows(y)-cols(x));

@ Ausdrucken @
format /rd 6,0; "Number of observations: "  rows(y);;
"  Degrees of freedom: " rows(y)-cols(x);
format /rd 6,3; "R2: " r2;; "  R2adj: " r2adj;

"     name  |  bp_hat   |    t-value | Corr. t-value | one-tail prob. | sig (two-tailed)";
?;
tvec =b./sqrt(diag(bvar));
tvec2=b./sqrt(diag(bvar2));
p1=cdfnc(abs(tvec2)); @ standard normal @

cvi=0;
do until cvi==rows(tvec);
cvi=cvi+1;
if cvi==1; let cv =" "; else; cv=cv|" "; endif;
endo;

cvi=0;
do until cvi==rows(tvec);
cvi=cvi+1;
if     p1[cvi] le (0.01)./2;                          cv[cvi]="***";
elseif p1[cvi] gt (0.01)./2 and p1[cvi] le (0.05)./2; cv[cvi]="**";
elseif p1[cvi] gt (0.05)./2 and p1[cvi] le (0.1)./2;  cv[cvi]="*";
endif;
endo;

omat=  seqa_
     ~ b
     ~ tvec
     ~ tvec2
     ~ p1
     ~ cv;

omat=sortc(omat,1);

mask=ones(1,5)~0;
let fmt[6,3] = "*.*lf"  8  0
               "*.*le" 16  3
               "*.*lf" 12  3
               "*.*lf" 12  3
               "*.*lf" 12  3
               "*.*s"   8  8;
call printfm(omat,mask,fmt);
if white.==1 and newest ne 1; "White(1980)-corrected standard errors"; endif;
if newest.==1; "Newey-West(1987)-corrected standard errors"; endif;

@ Wald-tests and F-tests: those with 1 on the diagonal will be eliminated @

"Wald- and F-test on significance of group of nonconstant regressors";
"Wald- and F-tests are based on corrected standard errors";
@ see Greene, 2nd ed., pages 188-189, 391, 392, 488, 491 @
"Test on all nonconstant regressors";
r=zeros(rows(b),rows(b));
@ those that are equal to one are kicked out @
r=diagrv(r,ones(rows(r),1));
r[rows(r),rows(r)]=0;
j=sumc(sumc(r)); @ number of restrictions @
gosub flag3;

"Wald- and F-test on significance of group of macro announcements";
@ see Greene, 2nd ed., pages 188-189, 391, 392, 488, 491 @
@ those that are equal to one are kicked out @
r=zeros(rows(b),rows(b));
rdiag=zeros(rows(b),1);
rdiag[1:12]=ones(12,1); @ to be dropped @
rdiag[14:rows(rdiag)-2]=ones(rows(rdiag)-14-1,1); @ to be dropped @
r=diagrv(r,rdiag);
j=sumc(sumc(r)); @ number of restrictions @
gosub flag3;

@ Jarque-Bera Test @

sig=sqrt(sige);
mu3=meanc(res.^3); mu4=meanc(res.^4);
lambd=((mu3.^2)./(6.*sig.^6)+((mu4-3.*sig.^4).^2)./(24.*sig.^8)).*rows(res);
p=cdfchic(lambd,2);
"Bera-Jarque: ";;
format /rd 8,3; lambd;; "  rhs-probability: " p;
if lambd ge 1e3; "for high BJ-values there may be numerical underflow"; endif;

endo; @ ig @

end;

flag3:
@ see Greene, page 392 @
i4=0; ex=zeros(rows(r),1);
do while i4<rows(r);
i4=i4+1;
if r[i4,.]==zeros(1,cols(r)); ex[i4]=1; endif;
endo;
r=delif(r,ex); @ Entferne Nullzeilen @
f=(1./j).*(r*b)'inv(r*bvar*r')*(r*b); @ F-test @
w=j.*f; @ Wald-test @

p1=cdfchic(w,j);
format /rd 5,0; "Wald-statistic: (";;j;;")";;
format /rd 8,3; w;; ",  rhs-probability:" p1;;
if p1 le 0.01; "***"; elseif p1 gt 0.01 and p1 le 0.05; "**";
elseif p1 gt 0.05 and p1 le 0.1; "*"; else; " "; endif;

v1=j; v2=rows(x)-rows(b);
p1=cdffc(f,v1,v2);
format /rd 5,0; "F-statistic: (";;v1;;";";;v2;;")";;
format /rd 8,3; f;; ",  rhs-probability:" p1;;
if p1 le 0.01; "***"; elseif p1 gt 0.01 and p1 le 0.05; "**";
elseif p1 gt 0.05 and p1 le 0.1; "*"; else; " "; endif;
return;



****************************************************************
****************************************************************
****************************************************************
****************************************************************
****************************************************************


(amemiya.prg)
@ Gauss program to perform econometric methods @
@ For a description of the program's functions see the program and comments below @
new; #lineson;
   
rao=1; @ 1: Rao test on heteroscedasticity @

save path=c:\data\strips\gauss;
output file=c:\prg\strips\amemiya.out reset;
str=datestr(0);
tr=timestr(0);

"GAUSS                            time  ";; tr;; "  date (US) ";; str;
?;"input file: amemiya.prg";

vnames="xx";
dtset1="c:\\data\\strips\\gauss\\dat2000";  @ data input for year 2000 @
dtset2="c:\\data\\strips\\gauss\\zinvzz";
dtset3="c:\\data\\strips\\gauss\\zinvzzzt";
dtset4="c:\\data\\strips\\gauss\\zinvzzz";

rob=0; @ 1: delete days on which a speech was given (not used in paper); 2: use core measures only @
match=0; @ shorten window for STRIPS/TIIS measure to make it comply with window of CM measure @

ig=0;
do while ig.<2;
ig=ig+1;

if match.==1 and ig.==2; break; endif; @ do STRIPS/TIIS measure only @

white=1;

vnames="xx";
dtset0="c:\\data\\strips\\gauss\\newey\\newey2";
dtset0="c:\\data\\yield\\gauss\\regress";
declare bvar2;

open f0=^dtset0;
dat=readr(f0,rowsf(f0));
f0=close(f0);
clear f0;

@ ----- define variables ----- @

if ig.==1;
  ?;?;"STRIPS (OTR) INFLATION MEASURE";?;
  if match./=1;
    dat=delif(dat,(dat[.,1].<=1).*(dat[.,2].<31).*(dat[.,3].==1997));
  elseif match.==1; @ use criterion used when ig.==2 @
    dat=delif(dat,dat[.,3].<1999);
    dat=delif(dat,(dat[.,1].<=1).*(dat[.,2].<5).*(dat[.,3].==1999));
  endif;
  y__=dat[.,4]; @ STRIPS inflation measure @
  /* format /rd 10,5; y__[1:10]; */
  dtset1="c:\\data\\strips\\gauss\\kernel\\res_otr";
/*
  dat[.,44]=missrv(dat[.,44],-9999);
  ex=dat[.,44]./=-9999;
  sumc(ex); wait;
*/
elseif ig.==2;
  ?;?;"CONSTANT MATURITY INFLATION MEASURE";?;?;
  dat=delif(dat,dat[.,3].<1999);
  dat=delif(dat,(dat[.,1].<=1).*(dat[.,2].<5).*(dat[.,3].==1999));
  y__=dat[.,5]; @ constant-maturity inflation measure @
  /* format /rd 10,5; y__[1:10]; */
  dtset1="c:\\data\\strips\\gauss\\kernel\\res_cm";
/*
  dat[.,44]=missrv(dat[.,44],-9999);
  ex=dat[.,44]./=-9999;
  sumc(ex); wait;
*/
endif;

@ compile dataset @

@ dat=dat[.,1:43 45 46]; @

@ use all variables bar the trading day indicator @
x__=dat[.,7:cols(dat)]~ones(rows(dat),1); @ standardized macro announcement variables~const @
let seqa_=1 2 3 6 5 8 9 7 10 11 12 13 33 35 34 15 17 16 18 25 19 20 21 23 24 26 29
          30 31 32 27 28 36 37 38 39 4 22 14 40 42;

@ cols(x__); rows(seqa_); wait; @

@ do not use the core numbers @
  x__=dat[.,7:12 14:35 37:39 41:cols(dat)]~ones(rows(dat),1);
                                     @ 12/13: CPI/Core; 35/36: PPI/Core;  39/40: Sales/ex.Auto @
  let seqa_=1 2 3 6 5 8   7 10 11 12 13 33 35 34 15 17 16 18 25 19 20 21 23 24 26 29
          30 31    27 28 36    38 39 4 22 14 40 42;
                                     @ 8/9: CPI/Core; 31/32: PPI/Core;  36/37: Sales/ex.Auto @

if rob.==2; @ use the core numbers @
  x__=dat[.,7:11 13:34 36:38 40:cols(dat)]~ones(rows(dat),1);
                                     @ 12/13: CPI/Core; 35/36: PPI/Core;  39/40: Sales/ex.Auto @
  let seqa_=1 2 3 6 5   9 7 10 11 12 13 33 35 34 15 17 16 18 25 19 20 21 23 24 26 29
          30    32 27 28    37 38 39 4 22 14 40 42;
                                     @ 8/9: CPI/Core; 31/32: PPI/Core;  36/37: Sales/ex.Auto @
endif;

x__=missrv(x__,0); @ turn empty cells into zeros @

y_=rev(y__);
x_=rev(x__);

@ obs are now stacked like this: 2003|2002|...|1997 @

@ now isolate missing observations (there was trading but we do not have the observation) @
@ note that days on which there is no trading have been eliminated from the dataset @

seqas=seqa(1,1,rows(dat)); @ index that includes missing observations @

y_=missrv(y_,-9);
ymiss=indexcat(y_,-9); @ positions of missing values @

y=delif(y_,y_.==-9);
x=delif(x_,y_.==-9);
seqam=delif(seqas,y_.==-9);

if rob.==1;
  y    =delif(y,x[.,cols(x)-1].==1);
  seqam=delif(seqam,x[.,cols(x)-1].==1);
  ymiss=ymiss|indexcat(x[.,cols(x)-1],1);
  x    =delif(x,x[.,cols(x)-1].==1);
endif;

wks=rows(seqas); @ number of obs in time, including missing values @
@ wks=rows(y); @
l=floor(4.*((wks./100).^(2./9))); @ lag for Ljung-Box and Newey-West @

/* ************** deal with speech/testimony indicator ************* */

z=        x[.,cols(x)-2:cols(x)]; @ do not exlude speech/testimony indicator: next to last column @
seqaz=seqa_[cols(x)-2:cols(x)];
zex_=(z[.,1].==0)+(z[.,2].==0);
zex=(zex_.==2);
@ z=delif(z,zex); @
@ create intercept indicator variable @
z=z[.,1 2]~zex~z[.,3]; @ two regressors, 0/1 intercept indicator, intercept @
seqaz=seqaz[1 2]|41|seqaz[3];
"number of nonzero observations:" rows(z)-sumc(z[.,3]);

x_=     x[.,1:cols(x)-2 cols(x)]; @ exlude speech/testimony indicator: next to last column @
seqax=seqa_[1:cols(x)-2 cols(x)];
@ create intercept indicator variable @
xex=(stdc(x_[.,1:cols(x_)-1]').==0);
x_=x_[.,1:cols(x_)-1]~xex~x_[.,cols(x_)];
seqax=seqax[1:rows(seqax)-1]|41|seqax[rows(seqax)];

/* ************************** core program ************************* */

y_=y;

@ modelling heteroscedasticity: includes speech/testimonry indicator @

z_=z;

@ dataset @

dst=y_~x_~z_;  @ z: regressors in the variance part of regression @
rx=cols(x_); @ number of regressors in main regression @
rz=cols(z_); @ number of regressors in variance regression @
nobs=rows(dst); @ number of observations @
clear y_,x_,z_;
/*
format /rd 8,6; dst;

omat=dst;
mask=ones(1,2);
let fmt[2,3] = "*.*lf"  8  3
               "*.*lf"  8  3;
call printfm(omat,mask,fmt);
  */
@ getting starting values @

y=dst[.,1]; @ dependent variable @
x=dst[.,2:(rx+1)]; @ main regression @
b=inv(x'x)*x'y;
res=y-x*b;
sige=(res'res)./rows(y);
bvar_=inv(x'x).*sige; @ uncorrected @

@ ------------------- White-correction ----------------- @

@ Greene, 2nd ed., page 391 @
hat=res.^2;
s0=zeros(cols(x),cols(x));
i=0;
do until i.==cols(x); @ across all regressors @
i=i+1;

s0[i,.]=(x[.,i]'.*hat')*x;

endo; @ i @
clear hat;

hat=diagrv(reshape(0,rows(x),rows(x)),res.^2);
/*
bvar=inv(x'x)*x'*hat*x*inv(x'x);
bvar=inv(x'x)*s0*inv(x'x);
*/
mat_=inv(x'x);
@ bvar=mat_*s0*mat_; White-corrected estimator of covariance matrix @

@ do loop for the matrix operation above @

mat_s0=reshape(0,rows(mat_),cols(s0));
i1=0;
do while i1<rows(mat_);
i1=i1+1;
i2=0;
do while i2<cols(s0);
i2=i2+1;
mat_s0[i1,i2]=mat_[i1,.]*s0[.,i2];
endo;
endo;

bvar=reshape(0,rows(mat_s0),cols(mat_));
i1=0;
do while i1<rows(mat_s0);
i1=i1+1;
i2=0;
do while i2<cols(mat_);
i2=i2+1;
bvar[i1,i2]=mat_s0[i1,.]*mat_[.,i2];
endo;
endo;

clear mat_,mat_s0;

tvec =b./sqrt(diag(bvar)); @ white @
tvec_=b./sqrt(diag(bvar_)); @ uncorrected @
p1=cdfnc(abs(tvec)); @ standard normal @

@---------------- end of White -----------------@

?; "OLS";
format /rd 6,0; "number of observations: "  rows(y);
"base regression";
"     coeff  |  s.e. |  t-value | White | one-tail prob. | sig (two-tailed)";
?;

cvi=0;
do until cvi==rows(tvec);
cvi=cvi+1;
if cvi==1; let cv =" "; else; cv=cv|" "; endif;
endo;

cvi=0;
do until cvi==rows(tvec);
cvi=cvi+1;
if     p1[cvi] le (0.01)./2;                          cv[cvi]="***";
elseif p1[cvi] gt (0.01)./2 and p1[cvi] le (0.05)./2; cv[cvi]="**";
elseif p1[cvi] gt (0.05)./2 and p1[cvi] le (0.1)./2;  cv[cvi]="*";
endif;
endo;

omat=  seqax
     ~ b
     ~ sqrt(diag(bvar))
     ~ tvec_ @ uncorrected @
     ~ tvec  @ white @
     ~ p1
     ~ cv;

omat=sortc(omat,1);

mask=1~ones(1,5)~0;
let fmt[7,3] = "*.*lf"  8  0
               "*.*le" 16  3
               "*.*le" 12  3
               "*.*lf" 12  3
               "*.*lf" 12  3               
               "*.*lf" 12  3
               "*.*s"   8  8;
call printfm(omat,mask,fmt);
format /le 14,3; "variance:" sige; @ wait; @

@---------------------- *************** ----------------------------@

@ Wald-tests and F-tests: those regressors that are removed from the model
  are assigned the value 1 on the first diagonal @
  
"Wald- and F-test on significance of group of nonconstant regressors";
@ see Greene, 2nd ed., pages 188-189, 391, 392, 488, 491 @
@ those that are equal to one are to be removed @
r=zeros(rows(b),rows(b));
r=diagrv(r,ones(rows(b),1));
r[rows(r),rows(r)]=0;
j=sumc(sumc(r)); @ number of restrictions @
call ftest(r,j,x,b,bvar);
/*
"Wald- and F-test on significance of group of year indicator variables";
@ see Greene, 2nd ed., pages 188-189, 391, 392, 488, 491 @
@ those that are equal to one are to be removed @
r=zeros(rows(b),rows(b));
rdiag=zeros(rows(b),1);
rdiag[cols(reg1~reg3)+1:cols(reg1~reg3~reg4)]=ones(cols(reg4),1); @ to be dropped @
r=diagrv(r,rdiag);
j=sumc(sumc(r)); @ number of restrictions @
call ftest(r,j,x,b,bvar);
endif;
*/
@ Jarque-Bera Test @

sig=sqrt(sige);
mu3=meanc(res.^3); mu4=meanc(res.^4);
lambd=((mu3.^2)./(6.*sig.^6)+((mu4-3.*sig.^4).^2)./(24.*sig.^8)).*rows(res);
p=cdfchic(lambd,2);
"Bera-Jarque: ";;
format /rd 8,3; lambd;; "  rhs-probability: " p;
if lambd ge 1e3; "for high BJ-values there may be numerical underflow"; endif;

@---------------------- analysis of residuals ------------------------------@

@ define variables--make sure constant regressor in in first spot, as required by Rao's score test) @


if rao.==1; @ Rao test on heteroscedasticity @

@ Rao's score test on homoscedasticity--Amemiya, 1985, pages 207-207, Equ. 6.5.33 @
/*
rao= (res.^2-sumc(res.^2)./rows(res))'
    *(z*inv(z'z)*z')
    *(res.^2-sumc(res.^2)./rows(res))
    ./(2.*(sumc(res.^2)./rows(res)).^2);
*/
 res2=     res.^2;
sres2=sumc(res.^2);
rres =rows(res);
invzz=inv(z'z);

@ --- create zinvzz ---- @

ih=0;
do while ih<rows(z);
ih=ih+1;
if     ih.==1; create f2=^dtset2 with ^vnames,cols(z),8;
               call writer(f2,z[ih,.]*invzz);
               f2=close(f2);
elseif ih.==2; open f2=^dtset2 for append;
               call writer(f2,z[ih,.]*invzz);
else;          call writer(f2,z[ih,.]*invzz);
endif;
endo;
f2=close(f2);
clear f2,invzz,ih;

@ --- create zinvzzz' ---- @

open f2=^dtset2 for read;
rw2=rowsf(f2);
zinvzz=readr(f2,rw2); @ z*invzz @
f2=close(f2);
clear f2;

ih=0;
do while ih<rows(z);
ih=ih+1;
if     ih.==1; create f3=^dtset3 with ^vnames,rows(z),8;
               call writer(f3,(zinvzz*z[ih,.]')');
               f3=close(f3);
elseif ih.==2; open f3=^dtset3 for append;
               call writer(f3,(zinvzz*z[ih,.]')');
else;          call writer(f3,(zinvzz*z[ih,.]')');
endif;
endo;
f3=close(f3);
clear f3,zinvzz,ih;

@ --- transpose zinvzzz' to zinvzzz --- @

open f3=^dtset3 for read;
rw3=rowsf(f3);
zinvzzzt=readr(f3,rw3); @ zinvzzz' @
f3=close(f3);
clear f3;

ih=0;
do while ih<rows(zinvzzzt);
ih=ih+1;
if     ih.==1; create f4=^dtset4 with ^vnames,cols(zinvzzzt),8;
               call writer(f4,zinvzzzt[.,ih]');
               f4=close(f4);
elseif ih.==2; open f4=^dtset4 for append;
               call writer(f4,zinvzzzt[.,ih]');
else;          call writer(f4,zinvzzzt[.,ih]');
endif;
endo;
f4=close(f4);
clear f4,zinvzzzt,ih;

@ --- transpose zinvzzz' to zinvzzz --- @

open f4=^dtset4 for read;
rw4=rowsf(f4);
zinvzzz=readr(f4,rw4); @ z*inv(z'z)*z' @
f4=close(f4);
clear f4;

sres2=sumc(res.^2);

rao= (res2-sres2./rres)'
    *zinvzzz
    *(res2-sres2./rres)
    ./(2.*(sres2./rres).^2);

format /ld 12,5; "Rao's score test (Amemiya, 1985, Eq. 6.5.33): " rao;; cdfchic(rao,cols(z)-1);

endif; @ rao.==1; @

@---------------------- residual regression ------------------------------@

alpha1=inv(z'z)*z'(res.^2); @ GQ estimator; Amemiya, 1985, Equ. 6.5.21: alpha1=inv(z'z)*z'(res.^2); @

@ covariance matrix for GQ estimator @

i=0;
do while i<rows(z);
i=i+1;
  if i.==1;
    c_=((z[i,.]*alpha1).^2)*(z[i,.]'*z[i,.]);
  else;
    c_=c_+((z[i,.]*alpha1).^2)*(z[i,.]'*z[i,.]);
  endif;
endo;

invzz=inv(z'z);
v1=2.*invzz*c_*invzz;
tvec=alpha1./sqrt(diag(v1));
p1=cdfnc(abs(tvec)); @ standard normal @

?; "residual regression"; 
"Goldfeld-Quandt estimators; Amemiya, 1985, Equ. 6.5.21: alpha1=inv(z'z)*z'(res.^2);";
"asympt. covariance matrix): Amemiya, 1985, Equ. 6.5.22";
"     coeff  |  s.e. |  t-value |one-tail prob. | sig (two-tailed)";
?;

cvi=0;
do until cvi==rows(tvec);
cvi=cvi+1;
if cvi==1; let cv =" "; else; cv=cv|" "; endif;
endo;

cvi=0;
do until cvi==rows(tvec);
cvi=cvi+1;
if     p1[cvi] le (0.01)./2;                          cv[cvi]="***";
elseif p1[cvi] gt (0.01)./2 and p1[cvi] le (0.05)./2; cv[cvi]="**";
elseif p1[cvi] gt (0.05)./2 and p1[cvi] le (0.1)./2;  cv[cvi]="*";
endif;
endo;

omat=  seqaz
     ~ alpha1
     ~ sqrt(diag(v1))
     ~ tvec
     ~ p1
     ~ cv;

omat=sortc(omat,1);

mask=1~ones(1,4)~0;
let fmt[6,3] = "*.*lf"   8  0
               "*.*le" 16  3
               "*.*le" 12  3
               "*.*lf" 12  3
               "*.*lf" 12  3               
               "*.*s"   8  8;
call printfm(omat,mask,fmt);

?;"Amemiya (1978) estimator; Amemiya, 1985, Equ. 6.5.30: alpha3=inv(z'*d1*z)*z'*d1*(res.^2);";
"asympt. covariance matrix: Amemiya, 1985, Equ. 6.5.31: avar=2.*inv(z'*d1*z);";

i=0;
do while i<rows(z);
i=i+1;
  if i.==1;
    c1=((z[i,.]*alpha1).^(-2))*(z[i,.]'*z[i,.]);
    c2=((z[i,.]*alpha1).^(-2))*(z[i,.]'*(res[i].^2));
    @ d1=((z[i,.]*alpha1).^(-2)); @
  else;
    c1=c1+((z[i,.]*alpha1).^(-2))*(z[i,.]'*z[i,.]);
    c2=c2+((z[i,.]*alpha1).^(-2))*(z[i,.]'*(res[i].^2));
    @ d1=d1|((z[i,.]*alpha1).^(-2)); @
  endif;
endo;

/* alternative calculation (tested)
d1=(z*alpha1).^-2; @ (column) vector @
d1=diagrv(eye(rows(d1)),d1);
alpha3=inv(z'd1*z)*z'd1*(res.^2); */

alpha3=inv(c1)*c2; @ Amemiya, 1985, Equ. 6.5.30: alpha3=inv(z'*d1*z)*z'*d1*(res.^2); @

avar=2.*inv(c1);  @ Amemiya, 1985, Equ. 6.5.31: avar=2.*inv(z'*d1*z); @

tvec=alpha3./sqrt(diag(avar));
p1=cdfnc(abs(tvec));

cvi=0;
do until cvi==rows(tvec);
cvi=cvi+1;
if cvi==1; let cv =" "; else; cv=cv|" "; endif;
endo;

cvi=0;
do until cvi.==rows(tvec);
cvi=cvi+1;
if     p1[cvi] le (0.01)./2;                          cv[cvi]="***";
elseif p1[cvi] gt (0.01)./2 and p1[cvi] le (0.05)./2; cv[cvi]="**";
elseif p1[cvi] gt (0.05)./2 and p1[cvi] le (0.1)./2;  cv[cvi]="*";
endif;
endo;

@ R-squared @
r2=1-(  sumc((res.^2-z*alpha3).^2)
       ./( (res.^2)'(res.^2)-rows(res.^2).*(meanc(res.^2).^2) ) );
r2adj=1-(1-r2) .* (rows(res.^2)-1) ./ (rows(res.^2)-cols(z));
format /rd 6,3; "R2: " r2;; "  R2adj: " r2adj;

omat= seqaz
     ~alpha3
     ~tvec
     ~p1
     ~cv;

omat=sortc(omat,1);

mask=1~ones(1,3)~0;
let fmt[5,3] = "*.*lf"  8  0
               "*.*le" 16  3
               "*.*lf" 12  3
               "*.*lf" 12  3
               "*.*s"   8  8;
call printfm(omat,mask,fmt);

endo; @ ig @

@------------------ procedure for Wald test and asymptotic F-Test ---------------------@

@ see Greene 2nd ed. page 392;  Greene 3rd ed. page 488 @
proc ftest(r,j,x,b,bvar);
  local ex,p1,f,v1,v2,i4,w;
  
  i4=0; ex=zeros(rows(r),1);
  do while i4<rows(r);
    i4=i4+1;
    if r[i4,.]==zeros(1,cols(r)); ex[i4]=1; endif;
  endo;

  r=delif(r,ex); @ Entferne Nullzeilen @
  f=(1./j).*(r*b)'inv(r*bvar*r')*(r*b); @ F-test @
  w=j.*f; @ Wald-test @
  /*
  p1=cdfchic(w,j);
  format /rd 5,0; "Wald-statistics: (";;j;;")";;
  format /rd 8,2; w;; ",  rhs-probability:" p1;;
  if p1 le 0.01; "***"; elseif p1 gt 0.01 and p1 le 0.05; "**";
  elseif p1 gt 0.05 and p1 le 0.1; "*"; else; " "; endif;
  */
  v1=j; v2=rows(x)-rows(b);
  p1=cdffc(f,v1,v2);
  format /rd 5,0; "F-statistics: (";;v1;;";";;v2;;")";;
  format /rd 8,5; f;; ",  rhs-probability:" p1;;
  
      if p1 le 0.01; "***"; elseif p1 gt 0.01 and p1 le 0.05; "**";
  elseif p1 gt 0.05 and p1 le 0.1; "*"; else; " ";
  endif;

retp(0);
endp;