This is the program file for Tables 2-7 in "The Asian Crisis and the Exposure of 
Large U.S. Firms" by William R. Emmons and Frank A Schmid in the 
January/February 2000 issue of the Federal Reserve Bank of St. Louis Review.
This file contains:  Gauss programs to create tables 2 through 7.

*********************
**Program Run Order**
*********************
a)  tab2.prg
b)  tab3_4.prg
c)  tab5.prg
d)  tab6.prg
e)  tab7.prg


@**************************************************************************************@
@Create Table 2@
@**************************************************************************************@
@ program tab2.prg @

new; #lineson;
save path=d:\data\asia1\gauss;
output file=c:\prg\asia4\tab2.out reset;
str=datestr(0);
tr=timestr(0);
"GAUSS                            time  ";; tr;; "  date (US) ";; str;
?;"input file: tab2.prg";

twodigit=1; @ 1: take two-digit codes; else: four-digit codes @

?;format /rd 4,0; screen on;

@ read data @
@ sort files by ticker symbol @

dataset="d:\\data\\asia1\\asc\\asia1";
open f1=^dataset;
rw=rowsf(f1);
asia1=readr(f1,rw); @ do not sort! (41 rows) @
f1=close(f1);

dataset="d:\\data\\asia1\\asc\\asia2";
open f1=^dataset;
rw=rowsf(f1);
asia2=readr(f1,rw); @ do not sort! (80 rows) @
f1=close(f1);

dataset="d:\\data\\asia1\\asc\\industry";
open f1=^dataset;
rw=rowsf(f1);
industry=sortcc(readr(f1,rw),2);
f1=close(f1);

if twodigit==1; industry[.,3]=floor(industry[.,3]./1e2); endif;

dataset="d:\\data\\asia1\\asc\\div97";
open f1=^dataset;
rw=rowsf(f1);
div97=sortcc(readr(f1,rw),2);
f1=close(f1);

dataset="d:\\data\\asia1\\asc\\div98";
open f1=^dataset;
rw=rowsf(f1);
div98=sortcc(readr(f1,rw),2);
f1=close(f1);

dataset="d:\\data\\asia1\\asc\\sales95";
open f1=^dataset;
rw=rowsf(f1);
sales95=sortcc(readr(f1,rw),1);
f1=close(f1);

dataset="d:\\data\\asia1\\asc\\sales96";
open f1=^dataset;
rw=rowsf(f1);
sales96=sortcc(readr(f1,rw),1);
f1=close(f1);

dataset="d:\\data\\asia1\\asc\\sales97";
open f1=^dataset;
rw=rowsf(f1);
sales97=sortcc(readr(f1,rw),1);
f1=close(f1);

dataset="d:\\data\\asia1\\asc\\sales98";
open f1=^dataset;
rw=rowsf(f1);
sales98=sortcc(readr(f1,rw),1);
f1=close(f1);

dataset="d:\\data\\asia1\\asc\\debt95";
open f1=^dataset;
rw=rowsf(f1);
debt95=sortcc(readr(f1,rw),1);
f1=close(f1);

dataset="d:\\data\\asia1\\asc\\debt96";
open f1=^dataset;
rw=rowsf(f1);
debt96=sortcc(readr(f1,rw),1);
f1=close(f1);

dataset="d:\\data\\asia1\\asc\\debt97";
open f1=^dataset;
rw=rowsf(f1);
debt97=sortcc(readr(f1,rw),1);
f1=close(f1);

dataset="d:\\data\\asia1\\asc\\debt98";
open f1=^dataset;
rw=rowsf(f1);
debt98=sortcc(readr(f1,rw),1);
f1=close(f1);

dataset="d:\\data\\asia1\\asc\\jpmweek";
open f1=^dataset;
rw=rowsf(f1);
jpm=readr(f1,rw);
f1=close(f1);

dataset="d:\\data\\asia1\\asc\\tb3week";
open f1=^dataset;
rw=rowsf(f1);
tb3=readr(f1,rw);
f1=close(f1);

dataset="d:\\data\\asia1\\asc\\ftweek";
open f1=^dataset;
rw=rowsf(f1);
ft=readr(f1,rw);
f1=close(f1);

@ create matrix of dummy variables  (second column: firm id) @

indus=unique(industry[.,3],1); @ SIC codes @
indu97=industry[.,2]~(industry[.,3].==indus');
indu98=indu97;
clear industry;

@ count # of SIC four-digit entries @

ex=(div98[.,3:cols(div98)]./=-9);
div98=div98[.,2]~(sumc(ex'));
ex=(div97[.,3:cols(div97)]./=-9);
div97=div97[.,2]~(sumc(ex'));

@ calculate sales and debt variables variables @

@ dummies @
s95=sales95;
s96=sales96;
s97=sales97;
d95=debt95;
d96=debt96;
d97=debt97;

@ fill in dummies @
k=0;
do while k<rows(sales98);
k=k+1;
if sales95[k,2] ne 12; @ 1995 @
       s95[k,3 5 7 9]= sales96[k,3 5 7 9].*(12-sales95[k,2])
                      +sales95[k,3 5 7 9].*    sales95[k,2];
           d95[k,2 3]=  debt96[k,2 3]    .*(12-sales95[k,2])
                      + debt95[k,2 3]    .*    sales95[k,2];
endif;
if sales96[k,2] ne 12; @ 1996 @
       s96[k,3 5 7 9]= sales97[k,3 5 7 9].*(12-sales96[k,2])
                      +sales96[k,3 5 7 9].*    sales96[k,2];
           d96[k,2 3]=  debt97[k,2 3]    .*(12-sales96[k,2])
                      + debt96[k,2 3]    .*    sales96[k,2];
endif;
if sales97[k,2] ne 12; @ 1997 @
       s97[k,3 5 7 9]= sales98[k,3 5 7 9].*(12-sales97[k,2])
                      +sales97[k,3 5 7 9].*    sales97[k,2];
           d97[k,2 3]=  debt98[k,2 3]    .*(12-sales97[k,2])
                      + debt97[k,2 3]    .*    sales97[k,2];
endif;
endo;
clear s95,s96,s97,d95,d96,d97;

@ reverse order (Newey-West) @

asia1=asia1[.,1 2]~(rev(asia1[.,3:cols(asia1)]'))';
asia2=asia2[.,1 2]~(rev(asia2[.,3:cols(asia2)]'))';
tb3  =             (rev(tb3')');
jpm  =             (rev(jpm')');
ft   =             (rev(ft')');

@ create variables @

ex=seqa(1,2,rows(sales95));
cap=sortcc(asia2[ex,.],2);
mtb=sortcc(asia2[ex+1,.],2);

@ check if files are of identical order @
@ a do-loop has to be used for characters @

i=0;
do while i<rows(sales95);
i=i+1;
crit1= cap[i,2];
if     mtb[i,2] $/= crit1
or  indu97[i,1] $/= crit1 or  indu98[i,1] $/= crit1
or   div97[i,1] $/= crit1 or   div98[i,1] $/= crit1
or sales95[i,1] $/= crit1 or sales96[i,1] $/= crit1
or sales97[i,1] $/= crit1 or sales98[i,1] $/= crit1
or  debt95[i,1] $/= crit1 or  debt96[i,1] $/= crit1
or  debt97[i,1] $/= crit1 or  debt98[i,1] $/= crit1;
"error in ordering, i:" i; "wait ";
wait; endif;
endo; clear crit1;

@ create vector of dmus @

dmuvec=seqa(1,1,rows(debt98));

@ eliminate Baker Hughes (BHI) because of imprecise sales data @
@ you may get a zero column in the industry dummy matrix @

ex=(sales95[.,1].$=="BHI"); @ look for ticker symbol @
/* ex=indexcat(ex,1); @ find line number @ */
 dmuvec=delif(dmuvec,ex);
    mtb=delif(mtb,ex);
    cap=delif(cap,ex);
 indu97=delif(indu97,ex);
 indu98=delif(indu98,ex);
  div97=delif(div97,ex);
  div98=delif(div98,ex);
sales95=delif(sales95,ex);
sales96=delif(sales96,ex);
sales97=delif(sales97,ex);
sales98=delif(sales98,ex);
 debt95=delif(debt95,ex);
 debt96=delif(debt96,ex);
 debt97=delif(debt97,ex);
 debt98=delif(debt98,ex);
@ rows(asia1)~rows(asia2); @
asia1=delif(asia1,asia1[.,2].$=="BHI");
asia2=delif(asia2,asia2[.,2].$=="BHI");
@ rows(asia1)~rows(asia2); @

@ eliminate firm id from indus and div files @

indu97=indu97[.,2:cols(indu97)];
indu98=indu98[.,2:cols(indu98)];
 div97= div97[.,2:cols(div97)];
 div98= div98[.,2:cols(div98)];

@ eliminate zero columns in indu97 and indu98 @

ex=(sumc(indu97|indu98)).<1;
if sumc(ex)>0;
indu97=delif(indu97',ex)';
indu98=delif(indu98',ex)';
endif;
clear ex;

@ define sales (lagged; fractions) @

s_usa = (sales97[.,3]./sales97[.,12])
       ~(sales96[.,3]./sales96[.,12]);
s_euro= (sales97[.,5]./sales97[.,12])
       ~(sales96[.,5]./sales96[.,12]);
s_asia= (sales97[.,7]./sales97[.,12])
       ~(sales96[.,7]./sales96[.,12]);
s_debt= ( debt97[.,2]./ debt97[.,3])
       ~( debt96[.,2]./ debt96[.,3]);
clear sales97,sales98,debt95,debt96,debt97,debt98;

@ check @

if maxc(maxc(s_debt))>1 or minc(minc(s_debt))<0; "error in s_debt -- wait";
wait; endif;
if maxc(maxc(s_asia))>1 or minc(minc(s_asia))<0; "error in s_asia -- wait";
wait; endif;

@ calculate herf of sales @

s_herf= sumc(((s_usa[.,1]~s_euro[.,1]~s_asia[.,1]).^2)')  @ 1997 @
       ~sumc(((s_usa[.,2]~s_euro[.,2]~s_asia[.,2]).^2)'); @ 1996 @

@ ranks @

rs_herf=s_herf; rcap=cap; rmtb=mtb; rs_usa=s_usa; rs_euro=s_euro;
rs_asia=s_asia;
rs_debt=s_debt; @ dummies @
k=2;
do while k<cols(cap);
k=k+1;
rcap[.,k]   =rankindx(cap[.,k],1);
rmtb[.,k]   =rankindx(mtb[.,k],1);
endo;
k=0;
do while k<cols(s_usa);
k=k+1;
rs_herf[.,k]=rankindx(s_herf[.,k],1);
rs_usa[.,k] =rankindx(s_usa[.,k],1);
rs_euro[.,k]=rankindx(s_euro[.,k],1);
rs_asia[.,k]=rankindx(s_asia[.,k],1);
rs_debt[.,k]=rankindx(s_debt[.,k],1);
endo;

@ alternative definition of ranks @

rs_usa = (sales96[.,3]./sales96[.,12])
        ~(sales95[.,3]./sales95[.,12]);
rs_euro= (sales96[.,5]./sales96[.,12])
        ~(sales95[.,5]./sales95[.,12]);
rs_asia= (sales96[.,7]./sales96[.,12])
        ~(sales95[.,7]./sales95[.,12]);
clear sales95,sales96;


@ calculate returns for individual firms, for S&P500 and for T-bill @
@ note: order has been reversed @

/* format /rd 6,2; tb3[.,1:2]; wait; */
rates=asia1[.,1 2]~ln(asia1[.,3:cols(asia1)-1]./asia1[.,4:(cols(asia1))]);
tb3 =ln( tb3[1:cols(tb3)-1] ./ tb3[2:cols(tb3)]);
rates=rates|(99~"tb3"~tb3);

@ calculate returns for indexes @
@ note: order has been reversed @

jpm_ =ln( jpm[1:cols(jpm)-1] ./ jpm[2:cols(jpm)]); @ weekly @
ft_  =ln(  ft[1,1:cols(ft)-1]  ./  ft[1,2:cols(ft)]);  @ weekly @

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

format /rd 6,0;
format /re 6,2;

div=div98~div97;
/*
i=0;
do while i<2;
i=i+1;

x=s_asia[.,i];

x=s_debt[.,i];
x=div[.,i];
/*
mtb_=mtb[.,3:(cols(mtb)-1)]';
x=mtb_[1:4,.]; x=reshape(x,cols(x)*rows(x),1);

cap_=cap[.,3:(cols(cap)-1)]';
x=cap_[5:8,.]; x=reshape(x,cols(x)*rows(x),1);
*/
gosub flag1;
endo;

end;

flag1:
format /rd 10,3;
minc(x);; median(x);; meanc(x);; maxc(x);; stdc(x);
return;
*/

format /rd 10,0; "div"; div;

end;


@**************************************************************************************@
@Tables 3 and 4@
@**************************************************************************************@
@ program Tab3_4.prg @

new; #lineson;
save path=d:\data\asia1\gauss;
output file=c:\prg\asia4\tab3_4.out reset;
str=datestr(0);
tr=timestr(0);
"GAUSS                            time  ";; tr;; "  date (US) ";; str;
?;"input file: tabs3_4.prg";

white0=0; @ 1: do white-correction @
numerair=0; @ 1: choose numeraire firm @
it=1e4; @ max. number of iterations @
contemp=0; @ 1: account for contemporaneous correlation @
twodigit=1; @ 1: take two-digit codes; else: four-digit codes @

?;format /rd 4,0; screen on;
if contemp==1; "contemporaneous correlation is accounted for";
else; "contemporaneous correlation is NOT accounted for"; endif;

sample=0; @ 1: take only a few firms @
save_=1; @ save res,x2,rx2,x104,rx104,b @
if sample==1; save_=0; endif; @ do not save @
if numerair==1; save_=0; endif;
iss1=0; @ must be equal to zero in any case @
iss2=0; @ 1: with jpm; 2: with ft @
do while iss2<1;
iss2=iss2+1;

iss=0; @ 1: test on significance in each year @
       @ 2: test on sign. difference across years @
do while iss<1;
iss =iss +1;
iss1=iss1+1;

@ convergence criterion @
if iss==1; crit=1e-10; elseif iss==2; crit=1e-10; endif;

@ read data @
@ sort files by ticker symbol @

dataset="d:\\data\\asia1\\asc\\asia1";
open f1=^dataset;
rw=rowsf(f1);
asia1=readr(f1,rw); @ do not sort! (41 rows) @
f1=close(f1);

dataset="d:\\data\\asia1\\asc\\asia2";
open f1=^dataset;
rw=rowsf(f1);
asia2=readr(f1,rw); @ do not sort! (80 rows) @
f1=close(f1);

dataset="d:\\data\\asia1\\asc\\industry";
open f1=^dataset;
rw=rowsf(f1);
industry=sortcc(readr(f1,rw),2);
f1=close(f1);

if twodigit==1; industry[.,3]=floor(industry[.,3]./1e2); endif;

dataset="d:\\data\\asia1\\asc\\div97";
open f1=^dataset;
rw=rowsf(f1);
div97=sortcc(readr(f1,rw),2);
f1=close(f1);

dataset="d:\\data\\asia1\\asc\\div98";
open f1=^dataset;
rw=rowsf(f1);
div98=sortcc(readr(f1,rw),2);
f1=close(f1);

dataset="d:\\data\\asia1\\asc\\sales95";
open f1=^dataset;
rw=rowsf(f1);
sales95=sortcc(readr(f1,rw),1);
f1=close(f1);

dataset="d:\\data\\asia1\\asc\\sales96";
open f1=^dataset;
rw=rowsf(f1);
sales96=sortcc(readr(f1,rw),1);
f1=close(f1);

dataset="d:\\data\\asia1\\asc\\sales97";
open f1=^dataset;
rw=rowsf(f1);
sales97=sortcc(readr(f1,rw),1);
f1=close(f1);

dataset="d:\\data\\asia1\\asc\\sales98";
open f1=^dataset;
rw=rowsf(f1);
sales98=sortcc(readr(f1,rw),1);
f1=close(f1);

dataset="d:\\data\\asia1\\asc\\debt95";
open f1=^dataset;
rw=rowsf(f1);
debt95=sortcc(readr(f1,rw),1);
f1=close(f1);

dataset="d:\\data\\asia1\\asc\\debt96";
open f1=^dataset;
rw=rowsf(f1);
debt96=sortcc(readr(f1,rw),1);
f1=close(f1);

dataset="d:\\data\\asia1\\asc\\debt97";
open f1=^dataset;
rw=rowsf(f1);
debt97=sortcc(readr(f1,rw),1);
f1=close(f1);

dataset="d:\\data\\asia1\\asc\\debt98";
open f1=^dataset;
rw=rowsf(f1);
debt98=sortcc(readr(f1,rw),1);
f1=close(f1);

dataset="d:\\data\\asia1\\asc\\jpmweek";
open f1=^dataset;
rw=rowsf(f1);
jpm=readr(f1,rw);
f1=close(f1);

dataset="d:\\data\\asia1\\asc\\tb3week";
open f1=^dataset;
rw=rowsf(f1);
tb3=readr(f1,rw);
f1=close(f1);

dataset="d:\\data\\asia1\\asc\\ftweek";
open f1=^dataset;
rw=rowsf(f1);
ft=readr(f1,rw);
f1=close(f1);

@ create matrix of dummy variables  (second column: firm id) @

indus=unique(industry[.,3],1); @ SIC codes @
indu97=industry[.,2]~(industry[.,3].==indus');
indu98=indu97;
clear industry;

@ count # of SIC four-digit entries @

ex=(div98[.,3:cols(div98)]./=-9);
div98=div98[.,2]~ln(sumc(ex'));
ex=(div97[.,3:cols(div97)]./=-9);
div97=div97[.,2]~ln(sumc(ex'));

@ calculate sales and debt variables variables @

@ dummies @
s95=sales95;
s96=sales96;
s97=sales97;
d95=debt95;
d96=debt96;
d97=debt97;

@ fill in dummies @
k=0;
do while k<rows(sales98);
k=k+1;
if sales95[k,2] ne 12; @ 1995 @
       s95[k,3 5 7 9]= sales96[k,3 5 7 9].*(12-sales95[k,2])
                      +sales95[k,3 5 7 9].*    sales95[k,2];
           d95[k,2 3]=  debt96[k,2 3]    .*(12-sales95[k,2])
                      + debt95[k,2 3]    .*    sales95[k,2];
endif;
if sales96[k,2] ne 12; @ 1996 @
       s96[k,3 5 7 9]= sales97[k,3 5 7 9].*(12-sales96[k,2])
                      +sales96[k,3 5 7 9].*    sales96[k,2];
           d96[k,2 3]=  debt97[k,2 3]    .*(12-sales96[k,2])
                      + debt96[k,2 3]    .*    sales96[k,2];
endif;
if sales97[k,2] ne 12; @ 1997 @
       s97[k,3 5 7 9]= sales98[k,3 5 7 9].*(12-sales97[k,2])
                      +sales97[k,3 5 7 9].*    sales97[k,2];
           d97[k,2 3]=  debt98[k,2 3]    .*(12-sales97[k,2])
                      + debt97[k,2 3]    .*    sales97[k,2];
endif;
endo;
clear s95,s96,s97,d95,d96,d97;

@ reverse order (Newey-West) @

asia1=asia1[.,1 2]~(rev(asia1[.,3:cols(asia1)]'))';
asia2=asia2[.,1 2]~(rev(asia2[.,3:cols(asia2)]'))';
tb3  =             (rev(tb3')');
jpm  =             (rev(jpm')');
ft   =             (rev(ft')');

@ create variables @

ex=seqa(1,2,rows(sales95));
cap=sortcc(asia2[ex,.],2);
mtb=sortcc(asia2[ex+1,.],2);

@ check if files are of identical order @
@ a do-loop has to be used for characters @

i=0;
do while i<rows(sales95);
i=i+1;
crit1= cap[i,2];
if     mtb[i,2] $/= crit1
or  indu97[i,1] $/= crit1 or  indu98[i,1] $/= crit1
or   div97[i,1] $/= crit1 or   div98[i,1] $/= crit1
or sales95[i,1] $/= crit1 or sales96[i,1] $/= crit1
or sales97[i,1] $/= crit1 or sales98[i,1] $/= crit1
or  debt95[i,1] $/= crit1 or  debt96[i,1] $/= crit1
or  debt97[i,1] $/= crit1 or  debt98[i,1] $/= crit1;
"error in ordering, i:" i; "wait ";
wait; endif;
endo; clear crit1;

@ create vector of dmus @

dmuvec=seqa(1,1,rows(debt98));

@ eliminate Baker Hughes (BHI) because of imprecise sales data @
@ you may get a zero column in the industry dummy matrix @

ex=(sales95[.,1].$=="BHI"); @ look for ticker symbol @
/* ex=indexcat(ex,1); @ find line number @ */
 dmuvec=delif(dmuvec,ex);
    mtb=delif(mtb,ex);
    cap=delif(cap,ex);
 indu97=delif(indu97,ex);
 indu98=delif(indu98,ex);
  div97=delif(div97,ex);
  div98=delif(div98,ex);
sales95=delif(sales95,ex);
sales96=delif(sales96,ex);
sales97=delif(sales97,ex);
sales98=delif(sales98,ex);
 debt95=delif(debt95,ex);
 debt96=delif(debt96,ex);
 debt97=delif(debt97,ex);
 debt98=delif(debt98,ex);
@ rows(asia1)~rows(asia2); @
asia1=delif(asia1,asia1[.,2].$=="BHI");
asia2=delif(asia2,asia2[.,2].$=="BHI");
@ rows(asia1)~rows(asia2); @

@ eliminate firm id from indus and div files @

indu97=indu97[.,2:cols(indu97)];
indu98=indu98[.,2:cols(indu98)];
 div97= div97[.,2:cols(div97)];
 div98= div98[.,2:cols(div98)];

@ eliminate zero columns in indu97 and indu98 @

ex=(sumc(indu97|indu98)).<1;
if sumc(ex)>0;
indu97=delif(indu97',ex)';
indu98=delif(indu98',ex)';
endif;
clear ex;

@ define sales (lagged; fractions) @

s_usa = (sales97[.,3]./sales97[.,12])
       ~(sales96[.,3]./sales96[.,12]);
s_euro= (sales97[.,5]./sales97[.,12])
       ~(sales96[.,5]./sales96[.,12]);
s_asia= (sales97[.,7]./sales97[.,12])
       ~(sales96[.,7]./sales96[.,12]);
s_debt= ( debt97[.,2]./ debt97[.,3])
       ~( debt96[.,2]./ debt96[.,3]);
/*
s_debt= ( debt97[.,3]./ (debt97[.,3]-debt97[.,2]))
       ~( debt96[.,3]./ (debt96[.,3]-debt96[.,2]));
s_debt=ln(s_debt);
*/
clear sales97,sales98,debt95,debt96,debt97,debt98;

@ check @

if minc(minc(s_debt))<0 or maxc(maxc(s_debt))>1; "error in s_debt -- wait";
wait; endif;
if maxc(maxc(s_asia))>1 or minc(minc(s_asia))<0; "error in s_asia -- wait";
wait; endif;

@ calculate herf of sales @

s_herf= sumc(((s_usa[.,1]~s_euro[.,1]~s_asia[.,1]).^2)')  @ 1997 @
       ~sumc(((s_usa[.,2]~s_euro[.,2]~s_asia[.,2]).^2)'); @ 1996 @

@ ranks period by period @

rs_herf=s_herf; rcap=cap; rmtb=mtb; rs_usa=s_usa; rs_euro=s_euro;
rs_asia=s_asia;
rs_debt=s_debt; @ dummies @
k=2;
do while k<cols(cap);
k=k+1;
rcap[.,k]   =rankindx(cap[.,k],1);
rmtb[.,k]   =rankindx(mtb[.,k],1);
endo;
k=0;
do while k<cols(s_usa);
k=k+1;
rs_herf[.,k]=rankindx(s_herf[.,k],1);
rs_usa[.,k] =rankindx(s_usa[.,k],1);
rs_euro[.,k]=rankindx(s_euro[.,k],1);
rs_asia[.,k]=rankindx(s_debt[.,k].*s_asia[.,k],1); @ !!! @
rs_debt[.,k]=rankindx(s_debt[.,k],1);
endo;

@ ranks (not period by period) @

vcap =reshape(cap,rows(cap).*cols(cap),1);
vmtb =reshape(mtb,rows(mtb).*cols(mtb),1);
vherf=reshape(s_herf,rows(s_herf).*cols(s_herf),1);
vusa =reshape(s_usa,rows(s_usa).*cols(s_usa),1);
veuro=reshape(s_euro,rows(s_euro).*cols(s_euro),1);
vasia=reshape(s_asia,rows(s_asia).*cols(s_asia),1);
vdebt=reshape(s_debt,rows(s_debt).*cols(s_debt),1);

rcap   =rankindx(vcap,1);
rmtb   =rankindx(vmtb,1);
rs_herf=rankindx(vherf,1);
rs_usa =rankindx(vusa,1);
rs_euro=rankindx(veuro,1);
rs_asia=rankindx(vasia.*vdebt,1); @ !!! @
rs_debt=rankindx(vdebt,1);

clear vcap,vmtb,vherf,vusa,veuro,vasia,vdebt;

rcap   =reshape(rcap,rows(cap),cols(cap));
rmtb   =reshape(rmtb,rows(mtb),cols(mtb));
rs_herf=reshape(rs_herf,rows(s_herf),cols(s_herf));
rs_usa =reshape(rs_usa,rows(s_usa),cols(s_usa));
rs_euro=reshape(rs_euro,rows(s_euro),cols(s_euro));
rs_asia=reshape(rs_asia,rows(s_asia),cols(s_asia));
rs_debt=reshape(rs_debt,rows(s_debt),cols(s_debt));

@ alternative definition of instruments: lags @

rs_usa = (sales96[.,3]./sales96[.,12])
        ~(sales95[.,3]./sales95[.,12]);
rs_euro= (sales96[.,5]./sales96[.,12])
        ~(sales95[.,5]./sales95[.,12]);
/*
rs_asia= (sales96[.,7]./sales96[.,12])
        ~(sales95[.,7]./sales95[.,12]); */
clear sales95,sales96;


@ calculate returns for individual firms, for S&P500 and for T-bill @
@ note: order has been reversed @

/* format /rd 6,2; tb3[.,1:2]; wait; */
rates=asia1[.,1 2]~ln(asia1[.,3:cols(asia1)-1]./asia1[.,4:(cols(asia1))]);
tb3_ =ln( tb3[1:cols(tb3)-1] ./ tb3[2:cols(tb3)]);
rates=rates|(99~"tb3"~tb3_);

@ calculate returns for indexes @
@ note: order has been reversed @

jpm_ =ln( jpm[1:cols(jpm)-1] ./ jpm[2:cols(jpm)])
          -rates[rows(rates),3:cols(rates)];
ft_  =ln(  ft[1,1:cols(ft)-1]  ./  ft[1,2:cols(ft)])
          -rates[rows(rates),3:cols(rates)];

    if iss2==1; air =  jpm_;
                air2=  jpm[1:cols(jpm)-1];
                air3=  jpm[1:cols(jpm)-1]./tb3[1:cols(tb3)-1];
elseif iss2==2; air =  ft_;
                air2=  ft[1,1:cols(ft)-1];
                air3=  ft[1,1:cols(ft)-1]./tb3[1:cols(tb3)-1];
endif;

air3=ln(air3./air3[cols(air3)]); @ normalize to unity & take log @

ex98=seqa(1,1,52);
ex97=ex98+52;
jpm97=meanc(jpm_[ex97]');
jpm98=meanc(jpm_[ex98]');
 ft97=meanc(ft_[ex97]');
 ft98=meanc(ft_[ex98]');
if save_==1; save jpm97,jpm98,ft97,ft98; endif;
format /le 10,3; (jpm97|jpm98)~(ft97|ft98);

@ # of dmus @

i=0; dmu=0;
if sample==1; firms=4; else; firms=40; endif;
do while i<firms; @ this number must read 40 ! for the full sample @
i=i+1;
ex=(dmuvec.==i); @ match firm id with index @
if sumc(ex) ge 1; dmu=dmu+1; endif;
endo;

if iss1>1; ?;
"-----------------------------------------------------------------------------";
?; endif;
if iss2==1; "jpm"; elseif iss2==2; "ft"; endif;
    if iss==1; "betas for <1998> and <1997>";
elseif iss==2; "betas for <1997 and 1998> difference of 1997 betas to them";
endif;

format /rd 4,0; "Number of equations: " dmu;

wks=cols(rates)-2; @ number of weekly observations @
l=floor(4.*((wks./100).^(2./9))); @ lag for Ljung-Box and Newey-West @

@ const has to be in last column @

i=0; ki=0;
if sample==1; firms=4; else; firms=40; endif;
do while i<firms; @ this number must read 40 ! for the full sample @
i=i+1;

ex=indexcat(dmuvec,i);
if sumc(ex) ge 1; ki=ki+1;

   @ blow up quarterly numbers @
   mtb_=ln((mtb[ex,3:(cols(mtb)-1)]').*.ones(13,1));
   cap_=ln((cap[ex,3:(cols(cap)-1)]').*.ones(13,1));
   rmtb_=  ((rmtb[ex,3:(cols(rmtb)-1)]').*.ones(13,1));  @ ranks @
   rcap_=  ((rcap[ex,3:(cols(rcap)-1)]').*.ones(13,1));  @ ranks @

   @ blow up annual numbers @
   s_herf_ = s_herf[ex,.]'.*.ones(52,1);
   s_usa_  =  s_usa[ex,.]'.*.ones(52,1);
   s_euro_ = s_euro[ex,.]'.*.ones(52,1);
   s_asia_ = s_asia[ex,.]'.*.ones(52,1);
   s_debt_ = s_debt[ex,.]'.*.ones(52,1);
   indust_ = (indu98[ex,.].*.ones(52,1))
            |(indu97[ex,.].*.ones(52,1));
      div_ = ( div98[ex,.].*.ones(52,1))
            |( div97[ex,.].*.ones(52,1));
   /* format /rd 6,2; i~s_asia_'; wait; */
   rs_herf_=rs_herf[ex,.]'.*.ones(52,1);
   rs_usa_ = rs_usa[ex,.]'.*.ones(52,1);
   rs_euro_=rs_euro[ex,.]'.*.ones(52,1);
   rs_asia_=rs_asia[ex,.]'.*.ones(52,1);
   rs_debt_=rs_debt[ex,.]'.*.ones(52,1);

   @ do not blow up quarterly numbers: get end of year numbers @
   mtbx=ln(mtb[ex,7 11]');
   capx=ln(cap[ex,7 11]');
   rmtbx=  (rmtb[ex,7 11]');  @ ranks @
   rcapx=  (rcap[ex,7 11]');  @ ranks @
   @ format /rd 10,4; exp(capx); wait; @

   @ do not blow up annual numbers: chose one number per year @
   s_herfx =  s_herf[ex,.];
   s_usax  =  s_usa[ex,.]; @ s_usax; wait; @
   s_eurox = s_euro[ex,.];
   s_asiax = s_asia[ex,.];
   s_debtx = s_debt[ex,.];
   industx =  indu98[ex,.]
            | indu97[ex,.];
      divx =   div98[ex,.]
            |  div97[ex,.];
   /* format /rd 6,2; i~s_asiax'; wait; */
   rs_herfx= rs_herf[ex,.];
   rs_usax = rs_usa[ex,.];
   rs_eurox=rs_euro[ex,.];
   rs_asiax=rs_asia[ex,.];
   rs_debtx=rs_debt[ex,.];

   @ form xx: one observation per year @

   @ 52 obs per year @
   x_  = s_herf_~ s_asia_~s_debt_~div_~mtb_~cap_~indust_;
   rx_ =rs_herf_~rs_asia_~s_debt_~div_~mtb_~cap_~indust_;
   @ 1 obs per year @
   xx_ = s_herfx'~ s_asiax'~s_debtx'~divx~mtbx~capx~industx;
   rxx_=rs_herfx'~rs_asiax'~s_debtx'~divx~mtbx~capx~industx;

   @ S&P 100 return minus risk-free rate @
   reg1_=( rates[rows(rates)-1,3:cols(rates)]
          -rates[rows(rates),3:cols(rates)])';

   @ for beta weighted by crisis variable @

   yvec1=zeros(rows(reg1_)./2,1)|ones(rows(reg1_)./2,1);  @ zero in 1998 @
   yvec2= ones(rows(reg1_)./2,1)|zeros(rows(reg1_)./2,1); @ zero in 1997 @
   if iss==1;
   reg1 = (reg1_.*yvec2)          @ 1998 @
         ~(reg1_.*yvec1)          @ 1997 @
         ~((reg1_.*air').*yvec2)  @ 1998 @
         ~((reg1_.*air').*yvec1); @ 1997 @
   elseif iss==2;
   reg1 = (reg1_@.*yvec2@)          @ 1997&1998 @
         ~(reg1_.*yvec1)          @ 1997 @
         ~((reg1_.*air')@.*yvec2@)  @ 1997&1998 @
         ~((reg1_.*air').*yvec1); @ 1997 @
   endif;

   @ for White test @
   reg1w= (reg1_.*yvec2)          @ 1998 @
         ~(reg1_.*yvec1)          @ 1997 @
         ~((reg1_.*air').*yvec2)  @ 1998 @
         ~((reg1_.*air').*yvec1); @ 1997 @
   rreg1=reg1;

   @ for graphs @
   g1_=meanc(s_asia_); g2_=meanc(s_debt_);

   z = mtb_~ cap_~(indust_.*reg1_.*air');
   rz=rmtb_~rcap_~(indust_.*reg1_.*air'); @ matrix for IV approach @
   clear reg1_;

   x = reg1;
   rx=rreg1; @ IV approach @

   @ White (no instruments enter the White test) @

       if iss==1; w_=reg1w;
   elseif iss==2; w_=reg1w@~z[.,1:2]@;
   endif;

   w=w_; clear reg1;
   l=0;
   do while l<cols(w_);
   l=l+1;
   k=0;
   do while k<l;
   k=k+1;
   w=w~(w_[.,l].*w_[.,k]);
   endo;
   endo;
   w=w~ones(rows(w),1)@~z[.,3:cols(z)]~(z[.,3:cols(z)].^2)@;

   @ firm return minus risk-free rate @
   y=(rates[ex,3:cols(rates)]-rates[rows(rates),3:cols(rates)])';

   if ki==1; xmat=     x;   ymat=     y; wmat=     w; zmat=     z;
             rzmat=rz;       rxmat=rx;
             dmat97=yvec1; dmat98=yvec2;
             g1=g1_; g2=g2_;
             x104=x_; rx104=rx_;
             x2 =xx_;  rx2 =rxx_;
             air_=air; air2_=air2; air3_=air3;
   else;     xmat=xmat|x;   ymat=ymat|y; wmat=wmat|w; zmat=zmat|z;
             rzmat=rzmat|rz; rxmat=rxmat|rx;
             dmat97=dmat97|yvec1; dmat98=dmat98|yvec2;
             g1=g1|g1_; g2=g2|g2_;
             x104=x104|x_; rx104=rx104|rx_;
             x2  =x2|xx_;  rx2  =rx2|rxx_;
             air_=air_|air; air2_=air2_|air2; air3_=air3_|air3;
   endif;
   clear x,y,w,z,rz,rx,x_,rx_,xx_,rxx_,d;
endif;
endo;

if iss==1;
if save_==1;
  if iss2==1;
    x2jpm=x2; rx2jpm=rx2; x104jpm=x104; rx104jpm=rx104;
    airjpm=air_; air2jpm=air2_; air3jpm=air3_;
    ymatjpm=ymat; xmatjpm=xmat; rxmatjpm=rxmat;
    save x2jpm,rx2jpm,x104jpm,rx104jpm,airjpm,air2jpm,air2jpm,air3jpm;
    save ymatjpm,xmatjpm,rxmatjpm;
  elseif iss2==2;
    x2ft =x2; rx2ft =rx2; x104ft =x104; rx104ft=rx104;
    airft=air_; air2ft=air2_; air3ft=air3_;
    ymatft=ymat; xmatft=xmat; rxmatft=rxmat;
    save x2ft,rx2ft,x104ft,rx104ft,airft,air2ft,air3ft;
    save ymatft,xmatft,rxmatft;
  endif;
endif;
endif;

if save_==1; save dmat97,dmat98; endif;

@ factors @

z=zmat; clear zmat;
rz=rzmat; clear rzmat;

names1="beta"|"beta2"|"beta3"|"beta4";
namesx=names1;
i=0;
do while i<dmu-1;
i=i+1;
namesx=namesx|names1;
endo;

names2="SIC";
i=0;
do while i<rows(indus)-1;
i=i+1;
names2=names2|"SIC";
endo;

/* ****************** assemble matrices ************************* */

i=0;
do while i<dmu;
i=i+1;
  x_=xmat[((i-1).*wks+1):(i.*wks),.]; @ one dmu @
      if i==1;   x=x_~zeros(rows(x_),((dmu-1).*cols(x_)));
  elseif i==dmu; x=x|(zeros(rows(x_),((dmu-1).*cols(x_)))~x_);
  else;          x=x|(zeros(rows(x_),((i-1).*cols(x_)))
                      ~x_
                      ~(zeros(rows(x_),((dmu-i).*cols(x_)))));
  endif;
  rx_=rxmat[((i-1).*wks+1):(i.*wks),.]; @ one dmu @
      if i==1;   rx=rx_~zeros(rows(rx_),((dmu-1).*cols(rx_)));
  elseif i==dmu; rx=rx|(zeros(rows(rx_),((dmu-1).*cols(rx_)))~rx_);
  else;          rx=rx|(zeros(rows(rx_),((i-1).*cols(rx_)))
                      ~rx_
                      ~(zeros(rows(rx_),((dmu-i).*cols(rx_)))));
  endif;
endo;
clear x_,rx_;

@ create numeraire firm (first firm serves as numeraire @
if numerair==1;
 x[.,1:4]= xmat;
rx[.,1:4]=rxmat;
endif;
@ clear xmat,rxmat; @

if iss<3;
  exogs=0;
  namesx=namesx;
else;
  rz=rz[.,3:cols(rz)]; z=z[.,3:cols(z)]; exogs=cols(z);
  namesx=@"LnMTB"|"LnCap"|"alpha"|@names2|namesx;
endif;

if iss>2;
  rx=rz~@ones(rows(x),1)~@rx; clear rz; @ ...with instruments @
   x= z~@ones(rows(x),1)~@ x; clear z;  @ add "factors", common alpha @
endif;

y=ymat; clear ymat;
w=wmat; clear wmat;

@ set up matrices for Wald-tests and F-tests @

if numerair==1;
  @ beta @
  rb=zeros(cols(x),cols(x));
  i=exogs+1;
  do while i<rows(rb)-3;
  i=i+4;
  rb[i,i]=1;
  endo;

  @ beta2 @
  rb2=zeros(cols(x),cols(x));
  i=exogs+2;
  do while i<rows(rb2)-2;
  i=i+4;
  rb2[i,i]=1;
  endo;

  @ beta3 @
  rb3=zeros(cols(x),cols(x));
  i=exogs+3;
  do while i<rows(rb3)-1;
  i=i+4;
  rb3[i,i]=1;
  endo;

  @ beta4 @
  rb4=zeros(cols(x),cols(x));
  i=exogs+4;
  do while i<rows(rb4);
  i=i+4;
  rb4[i,i]=1;
  endo;
else;
  @ beta @
  rb=zeros(cols(x),cols(x));
  i=exogs-3;
  do while i<rows(rb)-3;
  i=i+4;
  rb[i,i]=1;
  endo;

  @ beta2 @
  rb2=zeros(cols(x),cols(x));
  i=exogs-2;
  do while i<rows(rb2)-2;
  i=i+4;
  rb2[i,i]=1;
  endo;

  @ beta3 @
  rb3=zeros(cols(x),cols(x));
  i=exogs-1;
  do while i<rows(rb3)-1;
  i=i+4;
  rb3[i,i]=1;
  endo;

  @ beta4 @
  rb4=zeros(cols(x),cols(x));
  i=exogs;
  do while i<rows(rb4);
  i=i+4;
  rb4[i,i]=1;
  endo;
endif;

@ eliminate zero columns in x @

i=0; ex=zeros(cols(x),1);
do while i<cols(x);
i=i+1;
ex[i]=(x[.,i]==zeros(rows(x),1));
endo;
@ ex'; wait; @
x=(delif(x',ex))';
rx=(delif(rx',ex))';
namesx=delif(namesx,ex);
rb=(delif(rb',ex))';   rb=(delif(rb,ex));
rb2=(delif(rb2',ex))'; rb2=(delif(rb2,ex));
rb3=(delif(rb3',ex))'; rb3=(delif(rb3,ex));
rb4=(delif(rb4',ex))'; rb4=(delif(rb4,ex));

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

@ ohne kontemporaere Korrelation schaetzen @

b=inv(x'x)*x'y;
/* format /rd 6,4; b; wait; */
res=y-x*b;
resmat=reshape(res,dmu,wks); resmat=resmat';

sig=zeros(dmu,dmu); @ dummy-matrix; dmu: number of equations @
i=0;
do while i<dmu;
i=i+1;
j=0;
do while j<dmu;
j=j+1;
sig[i,j]=(resmat[.,i]'resmat[.,j])./rows(resmat);
endo;
endo;

imat=diagrv(zeros(wks,wks),ones(wks,1)); @ identity matrix @

if contemp==1; @ 1: account for cont. corr. @

   iti=0; break_=0; clear bold;
   do while iti<it; @ Iterationsloop @

   if break_==1; format /rd 2,0;
            "convergence achieved (";; iti;;")"; break; endif;
   iti=iti+1;

   b=inv(x'* (inv(sig).*.imat) *x)*x'* (inv(sig).*.imat) *y;
   bnew=b;
   if iti>1; critvec=abs((bnew-bold)./bold);
             format /rd 5,0;  "iteration # & iss" (iti-1)~iss;
             format /re 10,3; "   ...convergence:" sumc(critvec);
             if sumc(critvec)<crit; break_=1; endif;
   endif;
   bold=b;

   res=y-x*b;
   resmat=reshape(res,dmu,wks); resmat=resmat';

   sig=zeros(dmu,dmu); @ dummy-matrix; dmu: number of equations @
   i=0;
   do while i<dmu;
   i=i+1;
   j=0;
   do while j<dmu;
   j=j+1;
   sig[i,j]=(resmat[.,i]'resmat[.,j])./rows(resmat);
   endo;
   if iti==it; format /rd 2,0;
               "convergence not achieved (";; iti;;") -- wait"; endif;
   endo;

   endo; @ iti @

else; @ contemp==0 @
   sig=diagrv(zeros(dmu,dmu),diag(sig)); @ cont. corr. eliminated @
endif; @ contemp @

@ White test @

@ eliminate zero columns in w @

i=0; ex=zeros(cols(w),1);
do while i<cols(w);
i=i+1;
ex[i]=(w[.,i]==zeros(rows(w),1));
endo;
@ ex'; wait; @
w=(delif(w',ex))';

yw=res.^2;
bw=inv(w'w)*w'yw;
sst=sumc((yw-meanc(yw)).^2);
sse=(yw-w*bw)'(yw-w*bw);
r2=(sst-sse)./sst; clear yw,bw,sst,sse;
wstat=rows(res).*r2; clear r2;
pw1=cols(w)-1; clear w;
p1=cdfchic(wstat,pw1);
format /rd 5,0; "White-statistic: (";;pw1;;")";;
format /rd 8,3; wstat;; ",  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;

if p1<0.1; white=1; else; white=0; endif;

if white0.==1; white=0; endif;

wi=wstat~p1; clear wstat, p1;

if white==0;
bvar=inv(x'* (inv(sig).*.imat) *x); @ Greene, page 408, text below (17-8) @
endif;

@-------------------- Newey-West (1987) -----------------------@

@ Newey-West; see Greene, pp. 391, 422 @
@ Panel is balanced @
@ Obsverations are stacked like this: fonds1 1993
                                             1992
                                             ...
                                      fonds1 1993
                                             1992
                                             ...
                                         ... ...  @

@ Matrix s0 of White; see Greene, pp. 391, 422 @
@ SUR verlangt, dass statt der Matrix x, die Matrix xsurt und xsur
  verwendet wird @
xsur =   (inv(sig).*.imat)*x; @ [T        x equ] x [cols(x1) x equ] @
xsurt=x'*(inv(sig).*.imat);   @ [cols(x1) x equ] x [T        x equ] @

i6=0;
do while i6<rows(xsur); @ go over obs @
i6=i6+1;
if i6==1; s0=   (res[i6].^2).*(xsurt[.,i6]*xsur[i6.,]);
else;     s0=s0+(res[i6].^2).*(xsurt[.,i6]*xsur[i6.,]);
endif;
endo;
s0=s0./rows(x); @ divide by number of obs; Greene, page 391 @

@ --- newey-west (1987) --- @

@ generate a vector exvec that extracts the observations of the first fund @
exvec1=seqa(1,1,wks);

rj=zeros(l,2); @ fuer Ljung-Box; Greene, pages 426-27 @

i1=0; @ go through dmus (funds) @
do while i1<dmu;
i1=i1+1;

exvec=exvec1+(i1-1).*wks; @ choose a single fund @
/* exvec'; wait; */
res_=res[exvec];
regsur_ = xsur[exvec,.];
regsurt_=xsurt[.,exvec];

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]);

mod   =  w.*res_[i3].*res_[i3-i2]
       .*(regsurt_[.,i3]*regsur_[i3-i2,.]+regsurt_[.,i3-i2]*regsur_[i3,.]);
/* i2;; i3; wait; */
if i1==1 and i2==1 and i3==i2+1; s_star=mod; else; s_star=s_star+mod; endif;

endo; @ i3 @
endo; @ i2 @
endo; @ i1 @
s_star=s_star./rows(x); @ divide by number of obs @

@ 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);
p1=cdfchic(qprime,l);
format /rd 5,0; "Ljung-Box-statistic: (";;l;;")";;
format /re 12,3; qprime;; format /rd 8,4; ",  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;

if p1<0.1; serial=1; white=1; else; serial=0; endif; @ serielle Korrelation @

lj=qprime~p1; clear qprime, p1;

if white==1 and serial==1; white=0; serial=0; newest=1;
else; newest=0; endif;

if white ==1; s_ast=s0;        endif;
if serial==1; newest=1; endif; @ White-Teil muss mit verwendet werden! @
if newest==1; s_ast=s0+s_star; endif; @ Greene, page 422 @
clear s0,s_star;

@ Greene, 2nd ed., page 408 @
if serial==1 or white==1 or newest==1;
bvar=rows(x)
     .*inv(x'* (inv(sig).*.imat) *x) *s_ast *inv(x'* (inv(sig).*.imat) *x);
endif;

@-------------------- ENDE newey-west (1987) -----------------------@

@ format /rd 10,6; res'; wait; @
if iss==1;
if save_==1;
  if iss2==1;
    resjpm=res;
    save resjpm;
  elseif iss2==2;
    resft =res;
    save resft;
  endif;
endif;
endif;

@ R2 f�r jede Gleichung getrennt @

i=0; r2=seqa(1,1,dmu)~zeros(dmu,1);
do while i<dmu;
i=i+1;
y_=y[((i-1).*wks+1):(i.*wks)];
sig_=sig[i,i];
rrss=(y_-meanc(y_))'(y_-meanc(y_)); @ restricted @
urss=sig_.*wks; @ unrestricted @
r2[i,2]=1-urss./rrss; @ R-squared @
endo;

@ Ausdrucken @

format /rd 6,0; "Number of observations per equation: "  rows(y)./rows(sig);
if newest==1; "lag length L: " l; endif;
"# of equ. and R2";
omat=r2;
mask=ones(1,2);
let fmt[2,3] = "*.*lf"  6  0
               "*.*lf"  8  4;
call printfm(omat,mask,fmt);

?;"    variable    beta         t-value    one-tail probability";?;
tvec=b./sqrt(diag(bvar));
p1=cdfnc(abs(tvec));

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

cvi=0;
do while 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(1,1,rows(b))
     ~ namesx
     ~ b
     ~ tvec
     ~ cdfnc(abs(tvec))
     ~ cv;

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

if iss==1;
if save_==1;
  if iss2==1;
    bjpm=b~cv;
    save bjpm;
  elseif iss2==2;
    bft =b~cv;
    save bft;
  endif;
endif;
endif;

if white==1; ?;"Wald- and F-tests are based on White(1980)"; endif;
if newest==1; ?;"Wald- and F-tests are based on Newey-West(1987)"; endif;

@ Wald-tests and F-tests:
  1 on main diagonal for those regressors that are to be eliminated @

?;"Wald- and F-tests on significance of beta in system";
@ see Greene, 2nd ed., pages 188-189, 391, 392, 488, 491 @
r=rb;
j=sumc(sumc(r)); @ number of restrictions @
gosub flag3;

?;"Wald- and F-tests on significance of beta2 in system";
@ see Greene, 2nd ed., pages 188-189, 391, 392, 488, 491 @
r=rb2;
j=sumc(sumc(r)); @ number of restrictions @
gosub flag3;

?;"Wald- and F-tests on significance of beta3 in system";
@ see Greene, 2nd ed., pages 188-189, 391, 392, 488, 491 @
r=rb3;
j=sumc(sumc(r)); @ number of restrictions @
gosub flag3;

?;"Wald- and F-tests on significance of beta4 in system";
@ see Greene, 2nd ed., pages 188-189, 391, 392, 488, 491 @
r=rb4;
j=sumc(sumc(r)); @ number of restrictions @
gosub flag3;

flag4:?;

"Bera-Jarque test on normality:";

i=0; bj=seqa(1,1,dmu)~zeros(dmu,2);
do while i<dmu;
i=i+1;
res_=res[((i-1).*wks+1):(i.*wks)];
sig_=sqrt(sig[i,i]);
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);
bj[i,2 3]=lambd~p;
endo;

"# of equ., BJ-statistic and rhs-probability";
omat=bj;
mask=ones(1,3);
let fmt[3,3] = "*.*lf"  6  0
               "*.*le"  12 4
               "*.*lf"  8  4;
call printfm(omat,mask,fmt);

endo; @ iss @
endo; @ iss2 @

/* ******************************************************************** */

flag5: 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);
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 /re 10,3; w;; format /rd 5,3; ",  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 /re 10,3; f;; format /rd 5,3; ",  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;

@**************************************************************************************@
@Table 5@
@**************************************************************************************@
@ program tab5.prg @

new; #lineson;
save path=d:\data\asia4\gauss;
output file=c:\prg\asia4\tab5.out reset;
str=datestr(0);
tr=timestr(0);
"GAUSS                            time  ";; tr;; "  date (US) ";; str;
?;"input file: tab5.prg";

?;format /rd 4,0; screen on;

@ read data @

iss=0;
do while iss<2;
iss=iss+1;

if iss==1;

open f1=d:\data\asia1\gauss\bjpm.fmt;
rw=rowsf(f1);
b=readr(f1,rw);
f1=close(f1);

open f1=d:\data\asia1\gauss\x2jpm.fmt;
rw=rowsf(f1);
x=readr(f1,rw);
f1=close(f1);

open f1=d:\data\asia1\gauss\rx2jpm.fmt;
rw=rowsf(f1);
rx=readr(f1,rw);
f1=close(f1);

open f1=d:\data\asia1\gauss\jpm97.fmt;
rw=rowsf(f1);
w97=readr(f1,rw);
f1=close(f1);

open f1=d:\data\asia1\gauss\jpm98.fmt;
rw=rowsf(f1);
w98=readr(f1,rw);
f1=close(f1);

"jpm";

elseif iss==2;

open f1=d:\data\asia1\gauss\bft.fmt;
rw=rowsf(f1);
b=readr(f1,rw);
f1=close(f1);

open f1=d:\data\asia1\gauss\x2ft.fmt;
rw=rowsf(f1);
x=readr(f1,rw);
f1=close(f1);

open f1=d:\data\asia1\gauss\rx2ft.fmt;
rw=rowsf(f1);
rx=readr(f1,rw);
f1=close(f1);

open f1=d:\data\asia1\gauss\ft97.fmt;
rw=rowsf(f1);
w97=readr(f1,rw);
f1=close(f1);

open f1=d:\data\asia1\gauss\ft98.fmt;
rw=rowsf(f1);
w98=readr(f1,rw);
f1=close(f1);

?;?; "ft";

endif; @ iss @

ex98=seqa(3,4,39);
ex97=ex98+1;

iss1=0;
do while iss1<2;
iss1=iss1+1;

if iss1==1; ?; "1997";
db=b[ex97,1]~(w97.*b[ex97,1])~b[ex97,2];
elseif iss1==2; "1998";
db=b[ex98,1]~(w98.*b[ex98,1])~b[ex98,2];
endif;

omat=  seqa(1,1,rows(db))
     ~ db;

mask=ones(1,3)~0;
let fmt[4,3] = "*.*lf"  5  0
               "*.*le" 12  3
               "*.*le" 12  3
               "*.*s"   8  8;
screen on;
/*
call printfm(omat,mask,fmt);
*/
/* ************************** core program ************************* */

format /rd 6,0;
format /rd 6,3;

x=db[.,2]; gosub flag1;

endo; @ iss1 @
endo; @ iss @
end;

flag1:
format /le 10,2;
minc(x);; median(x);; meanc(x);; maxc(x);; stdc(x);
return;

@**************************************************************************************@
@Table 6@
@**************************************************************************************@
@ program tab6.prg @

/* instrumental variables approach */
/* pooled time series - cross section, balanced panel, Newey-West (1987);
   see Greene, 2nd ed. 1993, p. 379 (example) and p. 422) */

new; #lineson;
save path=d:\data\asia1\gauss;
output file=c:\prg\asia4\tab6.out reset;
str=datestr(0);
tr=timestr(0);
"GAUSS                            time  ";; tr;; "  date (US) ";; str;
?;"input file: tab6.prg";

?;format /rd 4,0; screen on;

@ read data @

iss=2;

if iss==1;

open f1=d:\data\asia1\gauss\bjpm.fmt;
rw=rowsf(f1);
y_=readr(f1,rw);
f1=close(f1);
/*
open f1=d:\data\asia1\gauss\ajpm.fmt;
rw=rowsf(f1);
y_=readr(f1,rw);
f1=close(f1);
y=y_[79:rows(y_),1];
*/
open f1=d:\data\asia1\gauss\x2jpm.fmt;
rw=rowsf(f1);
x=readr(f1,rw);
f1=close(f1);

open f1=d:\data\asia1\gauss\rx2jpm.fmt;
rw=rowsf(f1);
rx=readr(f1,rw);
f1=close(f1);
/*
"jpm";
*/
elseif iss==2;

open f1=d:\data\asia1\gauss\bft.fmt;
rw=rowsf(f1);
y_=readr(f1,rw);
f1=close(f1);
/*
open f1=d:\data\asia1\gauss\aft.fmt;
rw=rowsf(f1);
y_=readr(f1,rw);
f1=close(f1);
y=y_[79:rows(y_),1];
*/
open f1=d:\data\asia1\gauss\x2ft.fmt;
rw=rowsf(f1);
x=readr(f1,rw);
f1=close(f1);

open f1=d:\data\asia1\gauss\rx2ft.fmt;
rw=rowsf(f1);
rx=readr(f1,rw);
f1=close(f1);
/*
?;?; "ft";
*/
endif; @ iss @

ex=3|4; ex_=ex;
i=0;
do while i<(rows(y_)./4)-1;
i=i+1;
ex_=ex_|(ex+i.*4);
endo;
y=y_[ex_,1];

@ interaction between sales and debt @

 x[.,1:6]= x[.,1:6] .*x[.,3];
rx[.,1:6]=rx[.,1:6].*rx[.,3];
/*
 x= x.* x[.,3];
rx=rx.*rx[.,3];
*/
x= x[.,2 7:cols(x)-1];
rx=rx[.,2 7:cols(rx)-1];

names1=/*"s_herf"|
       */ "sales"
       /*|"debt"
       |"div"
       |"mtb"
       |"cap"*/;
/*
 x[.,cols(x)]= ones(rows(x),1);
rx[.,cols(x)]=ones(rows(rx),1);
*/
namesx=names1;

i=0;
do while i<18;
i=i+1;
namesx=namesx|"SIC";
endo;
/*
namesx=namesx|"const";
*/
wks=2; @ two periods @
dmu=rows(x)./wks;

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

@ lag length @
l=floor(4.*((wks./100).^(2./9))); @ lag for Ljung-Box and Newey-West @
if l==wks; l=l-1; endif;

@ OLS @

b=inv(rx'x)*rx'y;
res=y-x*b; @ residuals @
sige=(res'res)./rows(y); @ Varianz der Sch�tzung @
bvar=inv(rx'x)*(rx'rx)*inv(x'rx).*sige; @ Varianz der Sch�tzer b @

?; format /rd 6,0; "Number of observations: "  rows(y);

white=1;

@ Newey-West; see Greene, pp. 391, 422 @
@ Panel is balanced @
@ Obsverations are stacked like this: 1991 bank1
                                           bank2
                                           ...
                                      1990 bank1
                                           bank2
                                           ...
                                       ... ...  @

@ Matrix s0 of White; see Greene, pp. 391, 422 @
@ for s0, use instruments @
i1=0;
do while i1<rows(x);
i1=i1+1;
if i1==1; s0=   res[i1].^2.*(rx[i1,.]'*rx[i1,.]); @ xi': i-th row!!! @
else;     s0=s0+res[i1].^2.*(rx[i1,.]'*rx[i1,.]);
endif;
endo;
s0=s0./rows(x); @ dividiere durch Gesamtzahl der Beob. im Panel @

@ generate a vector exvec that extracts the observations of the first bank @
i2=0;
do while i2<wks; @ nobs @
i2=i2+1;
   if i2==1; exvec1=       (i2-1).*dmu+1;
   else;     exvec1=exvec1|((i2-1).*dmu+1); endif;
endo;

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

rj=zeros(l,2); @ fuer Ljung-Box; Greene, 2nd, pp. 426-7 @

i1=0; @ go through banks @
do while i1<dmu;
i1=i1+1;

exvec=exvec1+(i1-1); @ choose a single bank @
/* exvec'; wait; */
res_=res[exvec];
x_=rx[exvec,.]; @ use instruments @

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 i1==1 and i2==1 and i3==i2+1; s_star=mat; else; s_star=s_star+mat; endif;

endo; @ i3 @
endo; @ i2 @
endo; @ i1 @

@ calculate Ljung-Box @

rj[.,1]=rj[.,1]./sumc(res.^2); 
rj=(rj[.,1].^2)./rj[.,2]; 
qprime=sumc(rj).*(wks).*(wks+2);
p1lj=cdfchic(qprime,l);

if p1lj<0.05; serial=1; else; serial=0; endif; 

if serial==1; white=0; newest=1;
else; newest=0; endif;

if white ==1; s_ast=s0;        endif;
if newest==1; s_ast=s0+s_star; endif;
clear s0,s_star;

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

@-------------------- end of newey-west (1987) -----------------------@

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;
@ 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));
format /rd 6,3; "R2: " r2;; "  R2adj: " r2adj;
?;"     names         bp_hat       t-values  one-tail prob. sig (two-tailed)";
?;
tvec=b./sqrt(diag(bvar));
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;

omat=  seqa(1,1,rows(b))
     ~ namesx
     ~ b
     ~ tvec
     ~ p1
     ~ cv;

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

?;"If standard errors were corrected, Wald-tests and F-tests are based";
"on corrected standard errors";

?;"Wald- and F-test on all nonconstant regressors";
@ 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));
r=diagrv(r,ones(rows(b),1));
r[rows(r),rows(r)]=0;
j=sumc(sumc(r)); @ number of restrictions @
gosub flag3;

?;"Wald- and F-test on significance of industry dummies";
@ 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[2:cols(x)]=ones(18,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;

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

@**************************************************************************************@
@Table 7@
@**************************************************************************************@
@ program tab7.prg @

/* like newe1 but with bootstrapping */
/* instrumental variables approach */
/* pooled time series - cross section, balanced panel, Newey-West (1987);
   see Greene, 2nd ed. 1993, p. 379 (example) and p. 422) */

new; #lineson;
save path=d:\data\asia1\gauss;
output file=c:\prg\asia4\tab7.out reset;
str=datestr(0);
tr=timestr(0);
"GAUSS                            time  ";; tr;; "  date (US) ";; str;
?;"input file: tab7.prg";

?;format /rd 4,0; screen on;

runs=2.5e3; @2.5e3;  # of runs @
compute=0;

vnames="xx";
dtset1="d:\\data\\asia1\\gauss\\testvec1";

@ read data @

iss=1; @ 1: JPM; 2: FT/S&P @

if iss==1;

open f1=d:\data\asia1\gauss\bjpm.fmt;
rw=rowsf(f1);
y_=readr(f1,rw);
f1=close(f1);
/*
open f1=d:\data\asia1\gauss\ajpm.fmt;
rw=rowsf(f1);
y_=readr(f1,rw);
f1=close(f1);
y=y_[79:rows(y_),1];
*/
open f1=d:\data\asia1\gauss\x2jpm.fmt;
rw=rowsf(f1);
x=readr(f1,rw);
f1=close(f1);

open f1=d:\data\asia1\gauss\rx2jpm.fmt;
rw=rowsf(f1);
rx=readr(f1,rw);
f1=close(f1);
/*
"jpm";
*/
elseif iss==2;

open f1=d:\data\asia1\gauss\bft.fmt;
rw=rowsf(f1);
y_=readr(f1,rw);
f1=close(f1);
/*
open f1=d:\data\asia1\gauss\aft.fmt;
rw=rowsf(f1);
y_=readr(f1,rw);
f1=close(f1);
y=y_[79:rows(y_),1];
*/
open f1=d:\data\asia1\gauss\x2ft.fmt;
rw=rowsf(f1);
x=readr(f1,rw);
f1=close(f1);

open f1=d:\data\asia1\gauss\rx2ft.fmt;
rw=rowsf(f1);
rx=readr(f1,rw);
f1=close(f1);
/*
?;?; "ft";
*/
endif; @ iss @

ex=3|4; ex_=ex;
i=0;
do while i<(rows(y_)./4)-1;
i=i+1;
ex_=ex_|(ex+i.*4);
endo;
y=y_[ex_,1];

@ interaction between sales and debt @

 x[.,1:6]= x[.,1:6] .*x[.,3];
rx[.,1:6]=rx[.,1:6].*rx[.,3];
/*
 x= x.* x[.,3];
rx=rx.*rx[.,3];
*/
x= x[.,2 7:cols(x)-1];
rx=rx[.,2 7:cols(rx)-1];

names1=/*"s_herf"|
       */ "sales"
       /*|"debt"
       |"div"
       |"mtb"
       |"cap"*/;
/*
 x[.,cols(x)]= ones(rows(x),1);
rx[.,cols(x)]=ones(rows(rx),1);
*/
namesx=names1;

i=0;
do while i<18;
i=i+1;
namesx=namesx|"SIC";
endo;
/*
namesx=namesx|"const";
*/
wks=2; @ two periods @
dmu=rows(x)./wks;

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

tval1=inv(rx'x)*rx'y;

@ attach an index @

y =seqa(1,1,rows(y))~y;
x =seqa(1,1,rows(x))~x;
rx=seqa(1,1,rows(rx))~rx;

if compute==1;

@ uniform random variables @
rv=rndu(1e4,runs);

    if rows(x) gt 1e3 and rows(x) le 1e4; mult=1e4;
elseif rows(x) gt 1e2 and rows(x) le 1e3; mult=1e3;
elseif rows(x)                    le 1e2; mult=1e2;
else; "error(1)";
endif;

@ normalize random variables @
rv=trunc(rv.*mult);

@ start loop @

ex=-9.*ones(rows(x),runs); @ to be filled in @

i=0;
do while i<runs;
i=i+1;

j=0; j1=0;
do while j<rows(rv);
j=j+1;

ex_=(x[.,1].==rv[j,i]);
if sumc(ex_).==1; j1=j1+1; ex[j1,i]=rv[j,i];
if j1.==rows(x); break; endif;
endif;

endo;
endo;

if minc(minc(ex)) lt 0; "error(2) -- wait"; wait;
else; newe1aex=ex; save newe1aex; clear ex; endif;

endif; @ compute @

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

open f1=d:\data\asia1\gauss\newe1aex.fmt;
rw=rowsf(f1);
ex=readr(f1,rw); @ rows of x; # of cols = # or runs @
f1=close(f1);

i=0;
do while i<runs;
i=i+1;

subxh =seqa(2,1,cols(x)-1);
xh =submat(x,ex[.,i],subxh);
rxh=submat(rx,ex[.,i],subxh);
yh =submat(y,ex[.,i],2);

ex2=zeros(cols(xh),1);
j=0;
do while j<cols(xh);
j=j+1;
ex2[j]=(xh[.,j]==zeros(rows(xh),1));
endo;
if sumc(ex2).>0;
@ ex2'; wait; @
xh=(delif(xh',ex2))';
rxh=(delif(rxh',ex2))';
endif;

bh=inv(rxh'xh)*rxh'yh; @ IV approach @
bh=bh[1]; @ drop other coefficients because of zero columns,
              i.e, missing coefficients @

    if i==1; create f1=^dtset1 with ^vnames,cols(bh'),8;
             call writer(f1,bh');
             f3=close(f1);
elseif i==2; open f1=^dtset1 for append;
             call writer(f1,bh');
else;        call writer(f1,bh');
endif;

endo;

f1=close(f1);

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

open f1=^dtset1;
rw=rowsf(f1);
testvec1=sortc(readr(f1,rw),1);
f1=close(f1);

?; format /rd 4,0; "iss=" iss;

i=0;
do while i<cols(testvec1);
i=i+1;

@ confidence interval @

stdv=sqrt(sumc((testvec1[.,i]-meanc(testvec1[.,i])).^2)
          ./(rows(testvec1[.,i])-1));
format /re 10,3;
il1=1.6449; @ 90 percent @
il2=1.9600; @ 95 percent @
il3=2.576; @ 99 percent @

i1=0;
do while i1<3;
i1=i1+1;

    if i1.==1; il=il1; lbi=0.9;
elseif i1.==2; il=il2; lbi=0.95;
elseif i1.==3; il=il3; lbi=0.99;
endif;

(-il.*stdv)~tval1[i]~(il.*stdv);
@ wait; @

format /re 10,3;
lb=floor((1-lbi).*rows(testvec1[.,i]));
ub=floor(lbi.*rows(testvec1[.,i]));
testvec1=sortc(testvec1,i);
ci= (testvec1[lb,i]-meanc(testvec1[.,i]))~tval1[i]
   ~(testvec1[ub,i]-meanc(testvec1[.,i]));
ci;

@ wait; @
/*
format /re 10,3; testvec1'; @ "wait"; wait; @
hist(testvec1,11);
*/
endo; @ i1 @
endo; @ i @

end;