/*

This is the program file for Kevin L. Kliesen and Frank A. Schmid's 
   article "MACROECONOMIC NEWS AND REAL INTEREST RATES"
   in the Federal Reserve Bank of St. Louis's Review, March/April 2006 issue.
This file contains 3 Gauss programs: amemiya.prg, neweywest.prg and kernel01.prg

The programs use proprietary MMS data (not provided) and Fed Fund Surprise measures 
(available at: http://research.stlouisfed.org/publications/review/04/05/0405kkd.txt)

*/




//program part1 (amemiya.prg) This program is used to produce Table 3 "On the Run 10 year TIIS Yield and Data Surprises"

   
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=2; @ 1: delete days on which a speech was given (not used in paper); 2: use core measures only @

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

ln_=0; @ 1: take log first differences @
white=1;

vnames="xx";

dtset0="c:\\data\\strips\\gauss\\nabe"; @ Kuttner measure @
/* dtset0="c:\\data\\strips\\gauss\\nabe2"; @ Poole/Rasche measure @ */

declare bvar2;

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

@ compute the first differences @
@ obs are currently stacked like this: 1997|1998|...|2003 @

dat__=dat_[2:rows(dat_),.];

if ln_.==1;
  dat__[.,4 5]=ln(dat__[.,4 5]./dat_[1:rows(dat_)-1,4 5]);
else;
  dat__[.,4 5]=dat__[.,4 5]-dat_[1:rows(dat_)-1,4 5];
endif;

@ dat__[.,4 5]; wait; @

dat=dat__;
clear dat_,dat__;

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

y__4=dat[.,4]; @ change in 10-year TIIS yield @
y__5=dat[.,5]; @ change in 10-year TIIS yield @

if ig.==1;

    "OTR 10YR TIIS";
    y__=y__4;
    y_t=y__5; @ for match.==1 @
    

    dtset1="c:\\data\\strips\\gauss\\nwres1";

elseif ig.==2;

    "OTR 30YR TIIS";
     y__=y__5;
     y_t=y__4; @ for match.==1 @

    dtset1="c:\\data\\strips\\gauss\\nwres2";

elseif ig.==3;

    "OTR 30YR TIIS minus OTR 10YR TIIS";
     y__=y__5-y__4;
     y_t=y__4; @ for match.==1 @

    dtset1="c:\\data\\strips\\gauss\\nwres3";

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 8.1
          7 9 10 11 12 31 32 33 14 15 16 17 24 18 19 20 22 23 25 28 29
          30 30.1
          26 27
          34 34.1
          35 36 4 21 13
          41 42; @ intercept indicator; intercept @

@ 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 9 10 11 12 31 32 33 14 15 16 17 24 18 19 20 22 23 25 28 29
          30 
          26 27
          34 
          35 36 4 21 13
          41 42; @ intercept indicator; intercept @;

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
          8 
          7 9 10 11 12 31 32 33 14 15 16 17 24 18 19 20 22 23 25 28 29
          30 
          26 27
          34 
          35 36 4 21 13
          41 42; @ intercept indicator; intercept @
endif;

@ turn missing values into zero: for regressor matrix x only @

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

@ reverse order @

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

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

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); "number of nonzero observations"; rows(z); wait; */

@ create intercept indicator variable @
z=z[.,1 2]~zex~z[.,3]; @ two regressors (surprise, speech dummy), 0/1 intercept indicator, intercept @
z[.,1]=abs(z[.,1]); @ square the surprise variable or take absolute value @
@ format /rd 10,6; z; wait; @

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 @

"amemiya.prg finished gracefully";
end;

@------------------ 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;
//
//
//
////
//
//
//
//
//program part2 (neweywest.prg)
// This program is used for table 4 "Instrumental Variables Approach" and Table 5 "Uncertainty about Real Interest Rates"

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";

match=0; @ 1: use the same time window for 10- and 30-year yields @
ln_=0; @ 1: take log first differences @
_iv=0; @ 1: instrumental variables @
rob=2; @ 1: delete days on which a speech was given (not used in paper); 2: use core measures only @

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

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

vnames="xx";

dtset0="c:\\data\\strips\\gauss\\nabe"; @ Kuttner measure @
/* dtset0="c:\\data\\strips\\gauss\\nabe2"; @ Poole/Rasche measure @ */

declare bvar2, newest;

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

@ compute the first differences @
@ obs are currently stacked like this: 1997|1998|...|2003 @

dat__=dat_[2:rows(dat_),.];

if ln_.==1;
  dat__[.,4 5]=ln(dat__[.,4 5]./dat_[1:rows(dat_)-1,4 5]);
else;
  dat__[.,4 5]=dat__[.,4 5]-dat_[1:rows(dat_)-1,4 5];
endif;

@ dat__[.,4 5]; wait; @

dat=dat__;
clear dat_,dat__;

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

y__4=dat[.,4]; @ change in 10-year TIIS yield @
y__5=dat[.,5]; @ change in 10-year TIIS yield @

if ig.==1;

    "OTR 10YR TIIS";
    y__=y__4;
    y_t=y__5; @ for match.==1 @
    

    dtset1="c:\\data\\strips\\gauss\\nwres1";

elseif ig.==2;

    "OTR 30YR TIIS";
     y__=y__5;
     y_t=y__4; @ for match.==1 @

    dtset1="c:\\data\\strips\\gauss\\nwres2";

elseif ig.==3;

    "OTR 30YR TIIS minus OTR 10YR TIIS";
     y__=y__5-y__4;
     y_t=y__4; @ for match.==1 @

    dtset1="c:\\data\\strips\\gauss\\nwres3";

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 8.1
          7 9 10 11 12 31 32 33 14 15 16 17 24 18 19 20 22 23 25 28 29
          30 30.1
          26 27
          34 34.1
          35 36 4 21 13
          41 42; @ intercept indicator; intercept @

@ 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 9 10 11 12 31 32 33 14 15 16 17 24 18 19 20 22 23 25 28 29
          30 
          26 27
          34 
          35 36 4 21 13
          41 42; @ intercept indicator; intercept @;

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
          8 
          7 9 10 11 12 31 32 33 14 15 16 17 24 18 19 20 22 23 25 28 29
          30 
          26 27
          34 
          35 36 4 21 13
          41 42; @ intercept indicator; intercept @
endif;

@ turn missing values into zero: for regressor matrix x only @

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

@ reverse order @

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

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

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

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

y=delif(y_,y_.==minc(y_));
y_t=delif(y_t,y_.==minc(y_));
x=delif(x_,y_.==minc(y_));
seqam=delif(seqas,y_.==minc(y_));

if match.==1; @ do it for both dependent variables @
  y_t=missrv(y_t,minc(y_t)-1);
  y=delif(y,y_t.==minc(y_t));
  x=delif(x,y_t.==minc(y_t));
  seqam=delif(seqam,y_t.==minc(y_t));
endif;

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(seqam); @ 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; @

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

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

if _iv.==1;
  bvar=sige.*inv(z'x)*(z'*z)*inv(x'z);
else;
  bvar=inv(x'x).*sige; @ Varianz der Schaetzer b @
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;

if ismiss(ymiss) ne 1;
  resadd=zeros(rows(ymiss),1)~ymiss;
  xadd  =zeros(rows(ymiss),cols(x))~ymiss;
  res_=resnew|resadd;
  x_  =xnew|xadd;
else;
  res_=resnew;
  x_  =xnew;
endif;

if rows(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 @

"neweywest.prg finished gracefully";
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;
//
//
//
////
//
//
//
//
//prograp part 3 (kernel01.prg) 
// This program is used to produce Chart 1 "Distribution of Daily Changes in the Real Long Term Interest Rate"


@ May 20, 2004 @
/* To parts of this code, the followoing copyright applies:

**  (C) Copyright 1999 Gary King
**  All Rights Reserved.
**  http://GKing.Harvard.Edu, King@Harvard.Edu
**  Department of Government, Harvard University
*/
@------------------------------------------------------------------------@

@ Kernel estimation & optimization (if chosen) @
@ The kernel has to be Gaussian; select "0" for optimal bandwidth @
@ Optimal bandwidth is calculated as suggested by Silberberg (1986, eq. 3.28) @

@ new; #lineson; @
tstart=date;
library pgraph,optmum;
#include optmum.ext;

ig=1;
comp=0;

declare kernel01x,smth; @ for optimization (stores estimated parameters) @
vnames="xx";

dtset0="c:\\data\\strips\\gauss\\nabe"; @ Kuttner measure @
/* dtset0="c:\\data\\strips\\gauss\\nabe2"; @ Poole/Rasche measure @ */dtset1="c:\\data\\strips\\gauss\\wklyapril"; @ GSE data @

save path=c:\data\strips\gauss;

if comp.==1;
  output file=c:\prg\strips\kernel01.out reset;
endif;

"GAUSS                      time  ";; timestr(0);; "  date (US)";; datestr(0);
?;"input file: kernel01.prg";?;

if comp.==1;

dos del c:\data\strips\gauss\kernel01.*;

@ read the data @
match=1; @ 1: use the same time window for 10- and 30-year yields @
ln_=0; @ 1: take log first differences @
_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 @

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

@ compute the first differences @
@ obs are currently stacked like this: 1997|1998|...|2003 @

dat__=dat_[2:rows(dat_),.];

if ln_.==1;
  dat__[.,4 5]=ln(dat__[.,4 5]./dat_[1:rows(dat_)-1,4 5]);
else;
  dat__[.,4 5]=dat__[.,4 5]-dat_[1:rows(dat_)-1,4 5];
endif;

@ dat__[.,4 5]; wait; @

dat=dat__;
clear dat_,dat__;

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

if match.==1;

if ismiss(dat[.,4]).==1;
    dat[.,4]=missrv(dat[.,4],minc(dat[.,4])-1);
    dat=delif(dat,dat[.,4].==minc(dat[.,4]));
    endif;

if ismiss(dat[.,5]).==1;
    dat[.,5]=missrv(dat[.,5],minc(dat[.,5])-1);
    dat=delif(dat,dat[.,5].==minc(dat[.,5]));
    endif;

endif;

if ig.==1;

    "OTR 10YR TIIS";
    if ismiss(dat[.,4]).==1;
    dat[.,4]=missrv(dat[.,4],minc(dat[.,4])-1);
    dat=delif(dat,dat[.,4].==minc(dat[.,4]));
    endif;

    y__=dat[.,4]; @ change in 10-year TIIS yield @

    dtset1="c:\\data\\strips\\gauss\\nwres1";

elseif ig.==2;

    "OTR 30YR TIIS";
    if ismiss(dat[.,5]).==1;
    dat[.,5]=missrv(dat[.,5],minc(dat[.,5])-1);
    dat=delif(dat,dat[.,5].==minc(dat[.,5]));
    endif;

    y__=dat[.,5]; @ change in 10-year TIIS yield @

    dtset1="c:\\data\\strips\\gauss\\nwres2";

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 8
          7 9 10 11 12 31 32 33 14 15 16 17 24 18 19 20 22 23 25 28 29
          30 30
          26 27
          38 38
          39 40 4 21 13
          41 42; @ intercept indicator; intercept @

@ 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 9 10 11 12 31 32 33 14 15 16 17 24 18 19 20 22 23 25 28 29
          30 
          26 27
          38 
          39 40 4 21 13
          41 42; @ intercept indicator; intercept @;

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
          8 
          7 9 10 11 12 31 32 33 14 15 16 17 24 18 19 20 22 23 25 28 29
          30 
          26 27
          38 
          39 40 4 21 13
          41 42; @ intercept indicator; intercept @
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;

clear dat;

@ descriptive statistics @

  format /rd 6,4;
  @ minc(exp(y))~median(exp(y))~meanc(exp(y))~maxc(exp(y))~(stdc(exp(y)).*sqrt(rows(y)-1)./sqrt(rows(y))); @
  minc((y))~median((y))~meanc((y))~maxc((y))~(stdc((y)).*sqrt(rows(y)-1)./sqrt(rows(y)));

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

@ be sure to sort @

yk=sortc(y,1);
@ clear dat,y; @

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

@ calculate kurtosis and skewness @

sig=(meanc((yk-meanc(yk)).^2));
skew=(meanc((yk-meanc(yk)).^3))./(sig.^1.5);
skew_=cdfnc(abs(skew./sqrt(6./rows(yk))));

kurt=(meanc((yk-meanc(yk)).^4))./(sig.^2);
kurt_=cdfnc(abs((kurt-3)./sqrt(24./rows(yk))));

output on;
  format /rd 10,3; skew__=" "; kurt__=" ";
  if skew_ le 0.05; skew__="*"; endif; if skew_ le 0.025; skew__="**"; endif; if skew_ le 0.005; skew__="***"; endif;
  if kurt_ le 0.05; kurt__="*"; endif; if kurt_ le 0.025; kurt__="**"; endif; if kurt_ le 0.005; kurt__="***"; endif;
  "skewness (p-value):" skew;; skew_;; $skew__;
  "excess kurtosis (p-value):" kurt-3;; kurt_;; $kurt__;
output off;

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

@ kernel estimation @

{px1,py1}=dens(yk); 

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

cdfne=yk~seqa(0,1./(rows(yk)-1),rows(yk));

@-- procedure to compute objective function -- @

proc lpr(b); @ objective function @
     local f;
     
     f=sumc((cdfne[.,2]-cdfn((cdfne[.,1]-b[1])./b[2])).^2); @ cumulative normal distribution @ 
     
     retp(f);
endp;

@ starting values @

b0=meanc(yk)|stdc(yk);

optset;

_opgtol=reshape(1e-7,rows(b0),1); @ must be a vector; default: 1e-4 @
_opstmth = "NEWTON STEPBT"; 
__output=1;

call optprt(optmum_(&lpr,b0));

endif; @ comp.==1; @

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

open f1=c:\data\strips\gauss\kernel01x.fmt; @ estimated mean and standard deviation @
rw=rowsf(f1);
b=readr(f1,rw);
f1=close(f1);

b=b0; @ take the starting values instead @

graphset;
begwind;
makewind(9,6.855,0,0,0); @ Graph @

setwind(1);

grphstr="-w=20 -c=1 -cf=c:\\prg\\strips\\kernel01.eps";


_pltype={3 6 3};  @ line type @
_pcolor={8 8 8};  @ line color @

_pframe={0,0};
_pmcolor=8;
_paxht=.15; @ size characters of axes labels @
_ptitlht=.15; @ size of characters in title @
_pnotify=0; @ no screen activity while writing tkf file @
_plwidth={5}; @ _plwidth={1 3 1}; line width @
_pcsel=15;  @ white @
_pnum=2;    @ controls axes numbering @
_pnumht=0.15; @ size of axis numbering @
_plctrl={0 0 0};  @ frequency of symbols ; -1: symbols only @
_pstype=1;  @ 1: circle; symbol type @
_psymsiz=2; @ size of the symbols on the main curve @
_pltype=6; @ 6: solid @
fonts("simplex simgrma");

@ create box that covers up some of the axis info @

_pframe={0,0};

@ covers up part of the y-axis and the 0 @

makewind(0.94,1.13,0.1,0.1,0);

@ title("Return on Total Assets (1)"); @

ylabel("Probability Density");
xlabel("Logarithmic Excess Return");

  @ for graph @

  yk2=seqa(minc(yk),(maxc(yk)-minc(yk))./rows(yk),rows(yk));
  
  yk_=-((yk2-b[1])./b[2]).^2;
  py2_=(exp(yk_./2))./(b[2].*sqrt(2.*pi));

  @ for whiskers @
  
  yk_=-((yk-b[1])./b[2]).^2;
  py2=(exp(yk_./2))./(b[2].*sqrt(2.*pi));

scale(px1|yk,py1|py2);

if ig.==1;
  xtics(-0.15,0.25,0.05,2); ytics(-1,20,5,5); os6=0; os4=0;
elseif ig.==2;
  xtics(-0.25,0.25,0.05,2); ytics(-1,12,1,4); os6=0; os4=0;
endif;

os=ones(rows(yk),1);

  e1=unique(yk,1);
  os6=zeros(rows(yk),1);
  i2=0;
  do while i2<rows(e1);
    i2=i2+1;
    e2=indexcat(yk,e1[i2]);
    os6[e2]=sumc(e2).*ones(rows(e2),1);
  endo;
  os6=os6./(maxc(os6)./10)+os4;

_pline=   os
       ~ (os*6)              @ solid line @
       ~ yk                 @ x coordinates @
       ~ (os.*os4) @ beginning of whiskers; change with yticks @
       ~ yk
       ~ (os.*os6) @ ending of whiskers; change with ytics @
       ~  os        
      ~ (os*8)       @ color of line @
       ~(os*0);      @ thickness of line @

_plwidth={10};
xy(px1,py1); @ kernel @
_plwidth={5};
xy(yk2,py2_); @ normal pdf @

graphprt(grphstr);

endwind;

"kernel01.prg finished gracefully";
closeall;

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

/*  
**  (C) Copyright 1999 Gary King
**  All Rights Reserved.
**  http://GKing.Harvard.Edu, King@Harvard.Edu
**  Department of Government, Harvard University
**
**  Kernel density estimation
**
**   Format:    {px,py}=dens(y);
**
**   input:     y = Nx1 vector
**
**   output:    plot of y (horizontally) by its estd density f(y) (vertically).
**              Whiskers are also plotted along horizontal axis for each datum
**
**              px = x coordinates for the plot (_pts x 1 vector)
**              py = y coordinates for the plot (_pts x 1 vector)
**
**    Global values are as follows, and may be changed:
**   _smth=0;         smoothing parameter>0. the larger, the more smoothing.
**                    =0 for automatic calculation of density.
**                    =-1 for automatic calculation for description.
**
**   _strt=0;         x axis minimum. (if _strt=_endd, set at min and max).
**   _endd=0;         x axis maximum.
**   _pts=100;        number of points to plot
**   _whiskr=-1;      -1 = draw whiskers at 1/10th size;
**                    >0 = draw whiskers at size _whiskr;
**                     0 = do not draw whiskers
**                    Whiskers are always drawn at y==0, so you may have to
**                    use ytics to scale the y-axis.
**   _jitter=0;       0=nothing extra. #=add # jitter to each
**                      whisker if _whiskr/=0;.
**
**   let _kern=E;     kernel type. N=normal, E=Epanechnikov,
**                                 B=biweight, T=triangular, R=rectangular
**                                 TN=Doubly Truncated Normal
**  for option _kern=="TN",
**      _Tleft=_strt;   truncate distribution on left at _Tleft
**      _Tright=_endd;  truncate distribution on right at _Tright
**
**   _output = 1;     1 = print density plot (default); 0 = do not print.
**
**  Technical Reference:  
**    B.W. Silverman. 1986. _Density Estimation for Statistics
**    and Data Analysis_.  London: Chapman and Hall.
**
**  Application (code developed for):
**    King, Gary. "Constituency Service and Incumbency Advantage," British
**    Journal of Political Science, 21, 1 (January, 1991): 119-128, 
**    (replication dataset: ICPSR s1108)
*/
declare matrix _smth ?=0;
declare matrix _strt ?= 0;
declare matrix _endd ?= 0;
declare matrix _pts  ?= 2000;
declare matrix _Tleft ?= 0;
declare matrix _Tright ?= 0;
declare string _kern ?= "N";
declare matrix _whiskr ?= 0;
declare matrix _jitter ?= 0;
declare matrix _output ?= 1;
external proc code,xy,seqas_;
external matrix _pline;
@ Kernels @
@ arguments of kernels: z=input value, m=mean, h=smoothing parameter _smth @
proc kerneln(z,m,h);local res;z=(z-m)./h;      @ NORMAL kernel @
    res=(1./sqrt(2*pi))*exp(-(1/2).*(z^2));retp(res);endp;
proc kernelTN(z,m,h,Tleft,Tright);
    local t,tl,res,tr,zz;zz=(z-m)./h;          @ TRUNCATED NORMAL kernel @
    tl=(Tleft-m)./h; tr=(Tright-m)./h;
    t=((z.>Tleft).and (z.<Tright));
    res=t.*pdfn(zz)./(1-cdfn(tl)-(1-cdfn(tr)));
    retp(res);endp;
proc kernele(z,m,h);local a,res,t;z=(z-m)./h;  @ EPANECHNIKOV kernel @
    t=(abs(z).<sqrt(5)); a=code(t,sqrt(5)|1);
    res=t.*((3/4)*(1-(1/5).*(z^2))./a);retp(res);endp;
proc kernelb(z,m,h);local a,res,t;z=(z-m)./h;  @ BIWEIGHT kernel @
    t=(abs(z).<1);
    res=t.*((15/16)*((1-(z^2))^2));retp(res);endp;
proc kernelt(z,m,h);local a,res,t;z=(z-m)./h;  @ TRIANGULAR kernel @
    t=(abs(z).<1);
    res=t.*(1-abs(z));retp(res);endp;
proc kernelr(z,m,h);local a,res,t;z=(z-m)./h;  @ RECTANGULAR kernel @
    t=(abs(z).<1);
    res=t.*0.5;retp(res);endp;

@ creates two nx1 vectors to plot @
proc (2)=density(y,strt,endd,pts,h,kerna,Tleft,Tright);
    local i,t,px,py;    local kern:proc;
    kerna=upper(kerna);
    px=seqas_(strt,endd,pts);
    py=px;
    format /rdn 4,0;
    "Kernel Density Estimation";
    "Calculating ";;pts;;" Points: ";;
    i=1;do while i<=pts;
        if     (kerna$=="N"); t=kernelN(py[i,1],y,h)./h;
        elseif (kerna$=="E"); t=kernelE(py[i,1],y,h)./h;
        elseif (kerna$=="B"); t=kernelB(py[i,1],y,h)./h;
        elseif (kerna$=="T"); t=kernelT(py[i,1],y,h)./h;
        elseif (kerna$=="R"); t=kernelR(py[i,1],y,h)./h;
        elseif (kerna$=="TN");t=kernelTN(py[i,1],y,h,Tleft,Tright)./h;
        else; errorlog "_kernel specified incorrectly";stop;endif;
        py[i,1]=sumc( t )/rows(y);
        if i==10*floor(i/10);i;;endif;
    i=i+1;endo;?;retp(px,py);endp;

@ Kernel Density Estimate; sets defaults, calls density, then xy @
proc 2=dens(y);local strt,endd,h,kern,pts,px,py,smth,std,Tleft,Tright,os;
    if cols(y)/=1; errorlog "Argument must be a column vector";end;endif;
    strt=_strt;endd=_endd;
    if     strt>endd;errorlog "error: _strt>_endd";end;
    elseif strt==endd;strt=minc(y);endd=maxc(y);endif;
    pts=int(_pts);
    if pts<=2;errorlog "_pts must be greater than 2.  Try _pts=100;";endif;
    smth=_smth;
    if (smth<0).and(smth/=-1);errorlog "_smth must be -1 or > than 0";stop;
    elseif smth==0; py=sortc(y,1);
           std=minc((py[int(3*rows(py)/4)]-py[int(rows(py)/4)])/1.34|stdc(py));
           smth=0.9*std*(rows(py)^(-0.2));
           smth=(4./3).^(-0.2).*stdc(py).*(rows(py).^(-0.2)); 
    elseif smth==-1; py=sortc(y,1);
           std=minc((py[int(3*rows(py)/4)]-py[int(rows(py)/4)])/1.34|stdc(py));
           smth=0.9*std*(rows(py)^(-0.2))/2; endif;
    kern=_kern;
    Tleft=0;Tright=0;
    if kern$=="TN";
        if _Tleft==_Tright;  Tleft=strt;   Tright=endd;
        else;                Tleft=_Tleft; Tright=_Tright; endif;
    endif;
    {px,py}=density(y,strt,endd,pts,smth,kern,Tleft,Tright);
    std=rndu(rows(y),1)*_jitter-(_jitter/2);
    y=y+std;
    if _output==1;
    if _whiskr==-1;
        os=ones(rows(y),1);
        _pline=os~ (os*6)~ y~ (os*0)~ y~ (os*(maxc(py)/15))~ os~ (os*15)~(os*0);
        clear os;
    elseif _whiskr>0;
        os=ones(rows(y),1);
        _pline=os~ (os*6)~ y~ (os*0)~ y~ (os*_whiskr)~ os~ (os*15)~(os*0);
        clear os;
    endif;
        xy(px,py);
        format/rd 5,4;
        if rows(_smth)==1;?;"Smoothing Parameter: _smth=";;smth;endif;
    endif;
    retp(px,py);endp;
/*
**  (C) Copyright 1999 Gary King
**  All Rights Reserved.
**  http://GKing.Harvard.Edu, King@Harvard.Edu
**  Department of Government, Harvard University
**
**  creates vector of evenly spaced points 
**
**  y = seqas_(a,b,n);
**
**  y = vector of n points evenly spaced between a and b.  Y does not
**      include a or b
*/
proc seqas_(strt,endd,pts);
  local res;
  res=(endd-strt)/pts;
  res=seqa(strt+0.005*res,res,pts);
  retp(res);
endp;

/*
** optmum.src - General Nonlinear Optimization
** (C) Copyright 1988-2000 by Aptech Systems, Inc.
** All Rights Reserved.
**
** This Software Product is PROPRIETARY SOURCE CODE OF APTECH
** SYSTEMS, INC.    This File Header must accompany all files using
** any portion, in whole or in part, of this Source Code.   In
** addition, the right to create such files is strictly limited by
** Section 2.A. of the GAUSS Applications License Agreement
** accompanying this Software Product.
**
** If you wish to distribute any portion of the proprietary Source
** Code, in whole or in part, you must first obtain written
** permission from Aptech Systems.
**
**   Written by Ronald Schoenberg
**
**        CONTENTS                       LINE
**        --------                       ----
**        PROC OPTMUM                      29
**        Global Variables                 69
**        Using OPTMUM recursively        232
**        Source Code                     249
**
**-------------------**------------------**-------------------**-----------**
**-------------------**------------------**-------------------**-----------**
**
**   PROC OPTMUM
**
**   FORMAT
**          { x,f,g,retcode } = optmum_(&fct,x0)
**
**   INPUT
**
**        &fct - pointer to a procedure that computes the function to
**               be minimized.  This procedure must have one input
**               argument, a vector of parameter values, and one
**               output argument, the value of the function evaluated
**               at the input vector of parameter values.
**
**          x0 - vector of start values
**
**   OUTPUT
**          x - vector of parameters at minimum
**          f - function evaluated at x
**          g - gradient evaluated at x
**    retcode - return code:
**
**           0   normal convergence
**           1   forced exit
**           2   maximum number of iterations exceeded
**           3   function calculation failed
**           4   gradient calculation failed
**           5   Hessian calculation failed
**           6   step length calculation failed
**           7   function cannot be evaluated at initial parameter values
**           8   number of elements in the gradient vector inconsistent
**               with number of starting values
**           9   gradient function returned a column vector rather than
**               the required row vector
**          10   secant update failed
**          20   Hessian failed to invert
**
**
**-------------------**------------------**-------------------**-----------**
**-------------------**------------------**-------------------**-----------**
**
**   GLOBAL VARIABLES                                                 LINE
**
**   (default values in parantheses)
**  __title   -  string, title ("")                                    84
**  _opalgr   -  scalar, optimization algorithm (2)                    86
**  _opstep   -  scalar, selects type of step length (2)               93
**  _opshess  -  scalar or KxK matrix, selects starting Hessian (0)   105
**  _opfhess  -  KxK matrix, contains final Hessian                   111
**  _opmbkst  -  scalar, # of backsteps in computing steplength (10)  114
**  _opgtol   -  scalar, convergence tolerance for gradient (1e-5)    120
**  _opgdprc  -  scalar, pointer to gradient procedure (0)            124
**  _ophsprc  -  scalar, pointer to Hessian procedure (0)             143
**  _opgdmd   -  scalar, numerical gradient method (1)                158
**  _opparnm  -  Kx1 Char. vector, parameter names (0)                163
**  _opdfct   -  scalar, criterion for change in function (.001)      165
**  _opditer  -  scalar, # of iters to switch algorithms (20)         169
**  _opmiter  -  scalar, maximum number of iterations (1e+5)          173
**  _opmtime  -  scalar, maximum time in iterations in minutes (1e+5) 175
**  _oprteps  -  scalar, radius of random direction (0)               183
**  _opusrch  -  scalar, flag for user-controlled line search (0)     188
**  _opdelta  -  scalar, floor of Hessian eigenvalues in NEWTON (.1)  192
**  _opstmth  -  string, contains starting method ("")                196
**  _opmdmth  -  string, contains "middle" method ("")                204
**  _opkey    -  scalar, keyboard capture flag (1)                    210
**  _opgrdh   -  scalar, increment size for computing gradient (0)    219
**           
**   __title -  string, title of run
**
**    _opalgr -  scalar, indicator for optimization method:
**                = 1   SD (steepest descent - default)
**                = 2   BFGS (Broyden, Fletcher, Goldfarb, Shanno)
**                = 3   Scaled BFGS
**                = 4   Self-Scaling DFP (Davidon, Fletcher, Powell)
**                = 5   NEWTON (Newton-Raphson)
**                = 6   Polak-Ribiere Conjugate Gradient
**
**    _opstep - scalar, indicator determining the method for computing step
**              length.
**              = 1,  steplength = 1
**              = 2,  STEPBT  (default)
**              = 3,  golden steplength.
**              = 4,  Brent
**
**             Usually _opstep = 2 will be best.  If the optimization bogs
**             down try setting _opstep = 1 or 3.  _opstep = 3 will generate
**             slow iterations but faster convergence and _opstep = 1 will
**             generate fast iterations but slower convergence.
**
**   _opshess - scalar or KxK matrix, determines the starting hessian for
**              BFGS, DFP, and Newton methods.
**              = 0,  start with identity matrix  (default)
**              = 1,  computing starting hessian
**              = matrix,  user-defined starting hessian.
**
**   _opfhess - KxK matrix, contains the final Hessian, if one has been
**              calculated.  If the inversion of the Hessian fails during
**              NEWTON iterations this matrix can be analyzed for linear
**              dependencies which will suggest tactics for re-specifying
**              the model.
**
**  _opmbkst - scalar, maximum number of backsteps taken to find step length.
**             Default = 10.
**
**   _opgtol - scalar, convergence tolerance for gradient of estimated
**             coefficients.  Default = 1e-5.  When this criterion has been
**             satisifed OPTMUM will exit the iterations.
**
**  _opgdprc - scalar, pointer to a procedure that computes the gradient of the
**             function with respect to the parameters.  For example,
**             the instruction:
**
**                     _opgdprc=&gradproc
**
**             will tell OPTMUM that a gradient procedure exists as well
**             where to find it.  The user-provided procedure has a
**             single input argument, a vector of parameter values, and
**             a single output argument, a vector of gradients of the
**             function with respect to the parameters evaluated at the
**             vector of parameter values.  For example, suppose the
**             procedure is named gradproc and the function is a quadratic
**             function with one parameter: y=x^2+2*x+1, then
**
**                    proc gradproc(x); retp(2*x+2); endp;
**
**             Default = 0, i.e., no gradient procedure has been provided.
**
**  _ophsprc - scalar, pointer to a procedure that computes the hessian,
**             i.e., the matrix of second order partial derivatives of the
**             function with respect to the parameters.  For example, the
**             instruction:
**
**                    _ophsprc=&hessproc;
**
**             will tell OPTMUM that a procedure has been provided for the
**             computation of the hessian and where to find it.  The
**             procedure that is provided by the user must have a single
**             input argument, the vector of parameter values, and a single
**             output argument, the symmetric matrix of second order
**             derivatives of the function evaluated at the parameter
**             values.
**
**   _opgdmd - scalar, method for computing numerical gradient.
**                = 0, central difference
**                = 1, forward difference (default)
**                = 2, forward difference, Richardson Extrapolation
**
**  _opparnm - Kx1 character vector, parameter labels.
**
**   _opdfct - scalar, criterion for change in function which will cause
**             OPTMUM to switch algorithms when __design is nonzero.
**             Default = .001.
**
**  _opditer - scalar, criterion for maximum number of iterations before
**             switching algorithms when __design is nonzero.
**             Default = 20.
**
**  _opmiter - scalar, maximum number of iterations.  Default = 1e+5.
**
**  _opmtime - scalar, maximum time in iterations in minutes.
**              Default = 1e+5, about 10 weeks.
**
**  _oprteps - scalar, if _oprteps is set to a nonzero value (1e-2, say)
**             and all other line search methods fail then
**             OPTMUM will attempt a random direction with radius
**             determined by _oprteps.
**
**  _opusrch - scalar, if nonzero and if all other line search methods fail
**             OPTMUM will enter an interactive mode in which the
**             user can select a line search parameter.
**
**  _opdelta - scalar, used in NEWTON.  The eigenvalues of the Hessian
**             will be constrained to be greater than this value.
**             If set to zero the constraint will be disabled.
**
**  _opstmth - String, contains starting methods for algorithm, step
**             length, and Hessian.  For example,
**
**                  _opstmth = "NEWTON BRENT";
**
**             will cause OPTMUM to start with the NEWTON algorithm and
**             the BRENT step length method.
**
**  _opmdmth - String, contains "middle" methods for algorithm, step
**             length, and Hessian.  OPTMUM will switch to the methods
**             in _opmdmth either when the functions fails to change by
**             _opdfct, or when _opditer iterations have occured.
**
**    _opkey - scalar, flag controlling the keyboard capture feature.  During
**             the iterations OPTMUM will respond to the keyboard in order to
**             allow the user to modify global variables on the fly.  If
**             OPTMUM is being run recursively, however, (i.e., OPTMUM is
**             being called inside of the user procedure) then the keyboard
**             capture feature should be turned off for the version of
**             OPTMUM being called inside the user procedure which will
**             permit the outer version of OPTMUM to retain keyboard
**             control.
**
**   _opgrdh - scalar, increment size for computing gradient.
**
**
**-------------------**------------------**-------------------**-----------**
**-------------------**------------------**-------------------**-----------**
**
**   Calling OPTMUM recursively
**
**       The procedure provided by the user for computing the function to
**   to be minimized can itself call OPTMUM.  In fact the number of nested
**   levels is limited only by the amount of workspace memory.  Each level
**   will also contain its own set of global variables.  Thus nested versions
**   can have their own set of attributes and optimization methods.
**       It will be important to call OPTSET for all nested versions, and
**   generally if you wish the outer version of OPTMUM to retain control
**   over the keyboard you will need to set _opkey to zero for all the
**   nested versions.
**
*/

/*-------------------**------------------**-------------------**-----------**
**-------------------**------------------**-------------------**-----------*/

/*  SOURCE CODE  */

#include optmum.ext;
#include gauss.ext;

proc (4) = optmum_(fnct,x0);
    local x,f0,g,retcode;
    local Lopfhess,Lopitdta,LLoutput;

#ifUNIX
    LLoutput = __output /= 0;
#else
    LLoutput = __output;
#endif

    { x,f0,g,retcode,Lopfhess,Lopitdta } = _optmum(fnct,x0,
      _opalgr,
      _opdelta,
      _opdfct,
      _opditer,
      _opgdmd,
      _opgdprc,
      _opgrdh,
      _opgtol,
      _ophsprc,
      _opkey,
      _opmbkst,
      _opmdmth,
      _opmiter,
      _opmtime,
      _opmxtry,
      _opparnm,
      _oprteps,
      _opshess,
      _opstep,
      _opstmth,
      _opusrch,
      _opusrgd,
      _opusrhs,
      LLoutput,
      __title
      );

    _opfhess = Lopfhess;
    _opitdta = Lopitdta;

    retp(x,f0,g,retcode);
endp;



proc(0) = optset;
     gausset;
    _opalgr = 2;    /* optimization algorithm */
    _opparnm = 0;           /* parameter names */
    _opstep = 2;            /* selects type of step length */
    _opshess = 0;           /* selects starting hessian */
    _opmbkst = 10;          /* # of backsteps in computing steplength  */
    _opgtol = 1e-5;         /* convergence tolerance for gradient */
    _ophsprc = 0;           /* procedure to compute hessian */
    _opgdprc = 0;           /* procedure to compute gradient */
    _opgdmd = 0;            /* numerical gradient method */
    _opditer = 20;          /* # iters to switch algorithms for _opmdmth  */
    _opdfct = 0.001;        /* % change in function for _opmdmth */
    _opmiter = 1e+5;        /* maximum number of iterations */
    _opitdta = { 0,0,0 };
    _oprteps = .01;
    _opusrch = 0;
    _opdelta = .1;
    _opmxtry = 100;
    _opfhess = 0;   
    _opusrgd = 0;
    _opusrhs = 0;
    _opstmth = "";
    _opmdmth = "";
    _opkey = 1;
    _opgrdh = 0;

endp;

proc(4) = optprt(x,f,g,ret);
    local lbl,mask,fmt;
    print;
    call header("OPTMUM","",_op_ver);
    print;
    print "return code = " ftos(ret,"%*.*lf",4,0);
    if ret == 0;
       print "normal convergence";
    elseif ret == 1;
       print "forced termination";
    elseif ret == 2;
       print "maximum number of iterations exceeded";
    elseif ret == 3;
       print "function calculation failed";
    elseif ret == 4;
       print "gradient calculation failed";
    elseif ret == 5;
       print "Hessian calculation failed";
    elseif ret == 6;
       print "step length calculation failed";
    elseif ret == 7;
       print "function cannot be evaluated at initial parameter values";
    elseif ret == 8;
       print "number of elements in the gradient vector inconsistent";
       print "with number of starting values";
    elseif ret == 9;
       print "gradient function returned a column vector rather than";
       print "the required row vector";
    elseif ret == 11;
       print "maximum time exceeded";
    elseif ret == 20;
       print "Hessian failed to invert";
    endif;
    print;
    print "Value of objective function " ftos(f,"%*.*lf",15,6);
    print;
    print "Parameters    Estimates    Gradient";
    print "-----------------------------------------";

    if rows(g) /= rows(x);
         g = miss(zeros(rows(x),1),0);
    endif;

    if _opparnm $== 0;
        lbl = 0 $+ "P" $+ ftocv(seqa(1,1,rows(x)),2,0);
    else;
        lbl = _opparnm;
    endif;
    mask = 0~1~1;
    let fmt[3,3] = "-*.*s" 9 8 "*.*lf" 14 4 "*.*lf" 14 4;
    call printfm(lbl~x~g,mask,fmt);
    kernel01x=x; save kernel01x;    

    print;
    print "Number of iterations    " ftos(_opitdta[1],"%*.*lf",5,0);
    print "Minutes to convergence  " ftos(_opitdta[2],"%*.*lf",10,5);
    retp(x,f,g,ret);
endp;