clear;
echo off;
beep off;

% this is data from Fama and French
% ff data starts in 1926
% row 259 is 1948:1
% row 403  is 1960:1
% row 918 is 2002:12
% rows 750 is 1988:02
% row 846 is 1996:12
%load -ascii ~/matlab/cross/maxcv0.out;
%maxcv=maxcv0;
skips=0;
DEMEAN=2;
lags=6;
if skips==0;
%  cvfile='~/matlab/cross/maxcv0.out';
  cvfile='maxcv0.out';
maxcv=load(cvfile);
  else;
cvfile=['~/matlab/cross/maxcv0' int2str(skips*100) '.out'];
maxcv=load(cvfile);
end;



[dum1,dum2,dum3,dum4]=...
    textread('FFfac_m.asc','%d %f %f %f',910,'headerlines',1);
FF=[dum2 dum3 dum4];
factors=FF(403:906,:);
ffdate=dum1(403:906);


dir.data='swdata/';
%dir.main='~/matlab/riskdata/';
dir.main='';
mdfile='sw.inp';

dates.fyds = 1959;          % First Year Data Set
dates.fmds = 1;           % First Month Data Set
dates.lyds = 1996;          % Last Year Data Set
dates.lmds = 12;          % Last Month Data Set
dates.nfreq=12;

[series,tcode,rawdata,yy]=readdata(dir,mdfile,dates);

data1=trimr(yy,12,0);  % starts from 1960
dum1=1960+(1/12:1/12:47)';
[T,N]=size(data1);



% macro factors from DRI
load -ascii dridata/LHUR.dat;
ur=LHUR(13:456,3);
clear LHUR;
load -ascii dridata/IP.dat;
Mdate=IP(13:456,1)+IP(13:456,2)/12;
dip=mydiff(log(IP(13:456,3)),12);
clear IP;
load -ascii dridata/PCEND.txt;
dc=mydiff(log(PCEND(13:456,2)),12);
clear pcend;
load -ascii dridata/PUNEW.dat;
dp=mydiff(log(PUNEW(13:456,3)),12);
clear PUNEW;
load -ascii dridata/FYGM3.dat;
i3m=FYGM3(13:456,3);
clear FYGM3;
load -ascii dridata/FYGT10.dat;
i10y=FYGT10(13:456,3);
clear FYGT10;
term=i10y-i3m;
load -ascii dridata/FYBAAC.dat;
ibaa=FYBAAC(13:456,3);
clear FYBAAC;
risk=i10y-ibaa;
newfactors=lagn([dc dip dp term risk ur],1);
%newfactors=lagn(risk,1);

for subsample=1:5;

% from 1960 to 1996
if subsample==1;
  data=data1(1:444,:);
ffdate=dum1(1:444);
factor1=FF(1:444,:);
factor2=newfactors(1:444,:);
factor2=surprise(factor2,lags);
newdate=Mdate(1:444,:);
factors=[factor1 factor2];
delete('macro.6096');
diary('macro.6096');
end;



% 1960-1982
if subsample==2;
  data=data1(1:276,:);
ffdate=dum1(1:276);
factor1=FF(1:276,:);
factor2=newfactors(1:276,:);
factor2=surprise(factor2,lags);
newdate=Mdate(1:276,:);
factors=[factor1 factor2];
delete('macro.6082');
diary('macro.6082');
end;

% 83-96
if subsample==3;
  data=data1(277:444,:);
ffdate=dum1(277:444);
factor1=FF(277:444,:);
factor2=newfactors(277:444,:);
factor2=surprise(factor2,lags);
newdate=Mdate(277:444,:);
factors=[factor1 factor2];
delete('macro.8396');
diary('macro.8396');
end;

% 73-96
if subsample==4;
data=data1(157:444,:);
ffdate=dum1(157:444);
factor1=FF(157:444,:);
factor2=newfactors(157:444,:);
factor2=surprise(factor2,lags);
newdate=Mdate(157:444,:);
factors=[factor1 factor2];
delete('macro.7396');
diary('macro.7396');
end;


% 73-87
if subsample==5;
  data=data1(157:336,:);
ffdate=dum1(157:336);
factor1=FF(157:336,:);
factor2=newfactors(157:336,:);
factor2=surprise(factor2,lags);
newdate=Mdate(157:336,:);
factors=[factor1 factor2];
delete('macro.7387');
diary('macro.7387');
end;

x=[];
junk=[];
for i=2:cols(data);
  if sum(isnan(data(:,i)))>0;
        junk=[junk i];
    else;
    x=[x data(:,i)];
  end;
end;  


disp(sprintf('Demean %d',DEMEAN));
disp(sprintf('FF factors %d %d',[ffdate(1) ffdate(rows(ffdate))]));
disp(sprintf('macro variables %5f %5f',[newdate(1) newdate(rows(newdate),1)]));
[T,N]=size(x);
disp(sprintf('T= %d N = %d',T,N));
[maxcvL,maxcvH]=findcv(maxcv,.025,T);


rmax=10;
disp('detrmining number of factors');
for i=1:2;
  [ic1,chat1,fhat1]=nbpiid(x,rmax,i,DEMEAN);
  [ic2,chat2,fhat2]=nbplog(x,rmax,i,DEMEAN);
  disp([rmax ic1 ic2]);
end;  


[result,epshat,se_epshat,cr2,cr2se]=homotest(x,factors,6,maxcvH,skips,DEMEAN);
disp('homoskedasticity');
disp(result);
seup0=epshat+se_epshat;
sedn0=epshat-se_epshat;
%figure(1);
%subplot(2,1,1);
%plot(ffdate,[epshat(:,1) seup0(:,1) sedn0(:,1)]);
for i=1:cols(factors);
 disp(sprintf('%2d&%6.3f &%8.3f & %6.3f (%6.3f,%6.3f) & %6.3f & %6.3f(%6.3f,%6.3f) \\\\',...
    i,result.t(i),result.abs(i),result.r2(i),result.r2dn(i),result.r2up(i),...
      result.nir(i),cr2(i),cr2(i)-cr2se(i),cr2(i)+cr2se(i)));
  end;


[result,epshat,se_epshat]=heterotest(x,factors,6,maxcvH,skips,DEMEAN);
disp('heteroskedasticity');
disp(result)
seup1=epshat+se_epshat;
sedn1=epshat-se_epshat;
for i=1:cols(factors);
 disp(sprintf('%2d&%6.3f &%6.3f & %6.3f (%6.3f,%6.3f) & %6.3f & %6.3f(%6.3f,%6.3f)   \\\\',...
    i,result.t(i),result.abs(i),result.r2(i),result.r2dn(i),result.r2up(i),...
      result.nir(i),cr2(i),cr2(i)-cr2se(i),cr2(i)+cr2se(i)));

  end;
%figure(1);
%subplot(1,1,1);
%plot(ffdate, [epshat(:,5) seup1(:,5) sedn1(:,5)]);
%figure(2);
%subplot(1,1,1);
%plot(ffdate,[epshat(:,6) seup1(:,6) sedn1(:,6)]);
%figure(3);
%subplot(1,1,1);
%plot(ffdate,[epshat(:,9) seup1(:,9) sedn1(:,9)]);

%width0=seup0-sedn0;
%width1=seup1-sedn1;
% disp([i mean(x) max(x) min(x) std(x)]);


diary off;
end;
