% estimates favar by principal components
% confidence intervals are from a two-step bootstrap
% X is TxN, Y is TxM , F is TxK
% state equation:       [F(t)' Y(t)']'=B(L)[F(t-1)' Y(t-1)']'+v(t), v ~ N(0,Q), Q unrestricted
% X is Stock&Watson balanced panel less Y, observable factors
% Piotr Eliasz, 9/11/2002
clear all
close all
addpath('~/Dropbox/factor_models/BBE_Ddisk/BBE_Ddisk/Gibbs/')
addpath('~/Dropbox/factor_models/BBE_Ddisk/BBE_Ddisk/2step/')
addpath('~/Dropbox/factor_models/BBE_Ddisk/BBE_Ddisk')
% profile on -detail builtin

for num=1
    
    
    load US_data
    lab_english=textread('US_labels_english.m','%s','commentstyle','matlab','endofline','\n','delimiter','\n');
    rem_datasel=input('NEW SERIES? REMOVE PREVIOUS SERIES INFO? (0/1) ');
    if rem_datasel
        delete default_data_sel.mat
    end
    %% data selection: for subsample of countries/variables/periods
    if exist('default_data_sel.mat','file')==2;
        load default_data_sel.mat
    else
        datasel=1:length(onlynames);
        FFR='FFR';
        begsample=1;
        endsample=size(allFred,1);
    end
    %     datalabs=get_labels_excel('../dataUSFred/FXsub.xlsx','Monthly');
    sample=str2cell(sample);
    [datasel,FFR,savedef,begsample,endsample]=choose_subsample(lab_english,FFR,datasel,begsample,endsample,sample');
    if savedef
        save default_data_sel.mat datasel begsample endsample FFR -MAT
    end
    
    lab_english=lab_english(datasel,:);
    
    % data transformation
    %    -- Tcodes:
    %             1 Level
    %             2 First Difference
    %             3 Second Difference
    %             4 Log-Level
    %             5 Log-First-Difference
    %             6 Log-Second-Difference
    %             7 Detrend Log Using 1-sided HP detrending for Monthly data
    %             8 Detrend Log Using 1-sided HP detrending for Quarterly data
    %            16 Log-Second-Difference
    %            17 (1-L)(1-L^12)
    
    [tcode,name_var]=textread('US_labels.m','%u %s','commentstyle','matlab','endofline','\n','whitespace','\t');
    tcode=tcode(datasel);
    name_var=name_var(datasel,:);
    onlynames=onlynames(:,datasel);
    allFred=allFred(begsample:endsample,datasel);
    
    xdata=allFred;
    
    disp('variables used for US')
    disp([str2cell(num2str([1:length(name_var)]')),str2cell(num2str(tcode)),name_var,lab_english]);
    
    name_var_orig=name_var;
    name_var=lab_english;
    % Transform data to be approximately stationary
    for i_x = 1:size(xdata,2)   % Transform "X"
        
        if any(xdata(:,i_x)<0)&& tcode(i_x)==5
            error('Log of negative number in series # %i \n',i_x);
        end
        if any(isnan(xdata(:,i_x)))
            error('NaN in series # %i \n',i_x);
        end
        
        xtempraw(:,i_x) = transx(xdata(:,i_x),tcode(i_x)); %#ok<AGROW>
    end
    
    nsbalpanel=xtempraw(2:end,:);
    % need to standardize for PC:
    T=size(nsbalpanel,1);
    pos_FFR=find_incell(name_var,FFR);
    stdffr=std(nsbalpanel(:,pos_FFR));
    nsbalpanel=(nsbalpanel-repmat(mean(nsbalpanel,1),T,1))./repmat(std(nsbalpanel,1),T,1);
    if 0 % plot series to check stationarity
               zz=1;
        for kk=1:floor(size(nsbalpanel,2)/25)
            figure(kk);
            jj=1;
     
            while jj<=25
                subplot(5,5,jj)
                plot(nsbalpanel(:,zz))
                title(onlynames{zz})
                zz=zz+1;
                jj=jj+1;
            end
        end
        if size(nsbalpanel,2)-floor(size(nsbalpanel,2)/25)*25>0
            figure(kk+1);
            
            jj=1;
     
            while jj<=25
                subplot(5,5,jj)
                plot(nsbalpanel(:,zz))
                title(onlynames{zz})
                zz=zz+1;
                jj=jj+1;
            end
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%%% choose specification %%%%%%%%%%%%%%%%%%%%%%%%
    % Y=nsbalpanel(:,[16 107 77]); % this is IP, CPI and FFR
    Y=nsbalpanel(:,pos_FFR); % this is FFR
    xvars=setdiff([1:size(nsbalpanel,2)]',pos_FFR);
    X=nsbalpanel(:,xvars);%[nsbalpanel(:,1:pos_FFR-1),nsbalpanel(:,pos_FFR+1:end)];
    name_var=name_var(xvars,:);%[name_var(1:pos_FFR-1,:);name_var(pos_FFR+1:end,:)];
    % choose X from balanced panel (remember -- with Y excluded, so -1 after 77)
    xindex=[1:size(X,2)]';
    
    
    %varnames={'IP','CPI','3m TREASURY BILLS','5y TREASURY BONDS','M1','M2','MONETARY BASE','TOTAL INST RESERVES',...
    %        'NONBORROWED INST R','EXCHANGE RATE YEN','COMMODITY PRICE INDEX','CAPACITY UTIL RATE',...
    %        'PERSONAL CONSUMPTION','EMPLOYMENT','AVG WEEKLY HOURS',...
    %        'HOUSING STARTS','NEW ORDERS','NYSE COMPOSITE','AVG HOURLY EARNINGS',...
    %        'CONSUMER EXPECTATIONS'};
    
    
    slowindex=1;
    
    if exist('default_IRF.mat','file')==2;
        load default_IRF
    else
        checkedIRF=slowindex;
    end
    if max(slowindex)>length(name_var)
        slowindex=1;
        checkedIRF=slowindex;
    end
    [slowindex,numfac,lags,saveyes,checkedIRF,savedef]=slow_moving_list(name_var,slowindex,checkedIRF);
    if savedef
        save default_IRF.mat slowindex checkedIRF -MAT
    end
    % number of factors, lags in B(L):
    % numfac=[3];
    K=numfac(num);
    % lags=4;
    nrep1=10; % on observation equation
    nrep2=400; % on VAR equation
    nsteps=20;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    N=size(X,2);
    M=size(Y,2);
    
    % shock=[zeros(1,K+M-1) .078125]'; % in terms of standard deviation, identification is recursive
    shock=[zeros(1,K+M-1) .25/stdffr]'; % in terms of standard deviation, identification is recursive
    [FYresp,Xresp]=irfbootfac2step(Y,lags,X,xindex,slowindex,K,shock,nsteps-1,nrep1,nrep2);	% nsteps x equations x error bands
    FYresp=permute(FYresp,[2,1,3]);		% equations x nsteps x error bands
    Xresp=permute(Xresp,[2,1,3]);		% equations x nsteps x error bands
    
    VARresp=[squeeze(Xresp(1,:,1))' squeeze(Xresp(2,:,1))' squeeze(FYresp(end,:,1))'];
    
    name=['FAVARresp' int2str(K)];
    
    eval(['save ' name '.mat VARresp']);         % Saving results to output file: IRF's,
    %
    %
    % names={'FY','X'};
    % for n = 1:2
    %     save([mydir sprintf('%srespni_%d_%d_%d_%d.mat',names{n},K,lags,nrep1,nrep2)],...
    %         sprintf('%sresp',names{n}),'-mat');
    % end
    %clear  % when computing for many K's at one execution
end
% profile report


% Plotting Impulse Response Functions
% disp('CHOOSE WHETHER TO PRINT IN THE POP-UP MENU')
% saveyes=menu('SAVE IRFs','YES','NO');
resp=Xresp;

% index=[16 5;107 5;77 1;80 1;91 5;95 5;96 5;97 5;74 5;101 1;17 1;49 5;32 5;46 1;54 4;62 1;66 5;118 5;119 1];
% varnames={'IP','CPI','3m TREASURY BILLS','5y TREASURY BONDS','M1','MONETARY BASE','TOTAL INST RESERVES',...
%         'NONBORROWED INST R','EXCHANGE RATE YEN','COMMODITY PRICE INDEX','CAPACITY UTIL RATE',...
%         'PERSONAL CONSUMPTION','EMPLOYMENT','AVG WEEKLY HOURS',...
%         'HOUSING STARTS','NEW ORDERS','NYSE COMPOSITE','AVG HOURLY EARNINGS',...
%         'CONSUMER EXPECTATIONS'};


% index=[xindex,[tcode(1:pos_FFR-1);tcode(pos_FFR+1:end)]];
index=checkedIRF;

varnames=regexprep(name_var(checkedIRF,:),'\_','\\_');%([1:pos_FFR-1,pos_FFR+1:end],1);

nsteps=20;
n=size(resp,1);

% % transform back to levels
% for i=1:n
%     if index(i,2)==4
%         resp(i,:,:)=exp(resp(i,:,:))-ones(1,nsteps,3);
%     elseif index(i,2)==5
%         resp(i,:,:)=exp(cumsum(resp(i,:,:),2))-ones(1,nsteps,3);
%     end
% end
resp_old=resp;
resp=resp(checkedIRF,:,:);
s=1:nsteps; z=zeros(nsteps,1);	% to plot level zero
zz=0;
numser=size(resp,1);
fullfig=floor(numser/(15));
resid=numser-fullfig*15;
for jj=1:fullfig+min(1,resid);
    fg(jj)=figure;
    
    subplot(4,4,1)
    plot(s,[z squeeze(FYresp(end,:,1))'],'LineWidth',1.5);
    hold on
    h=ciplot(squeeze(FYresp(end,:,2))',squeeze(FYresp(end,:,5))',s,'r');
    set(h,'EdgeAlpha',.00)
    set(h,'FaceAlpha',.1)
    h=ciplot(squeeze(FYresp(end,:,3))',squeeze(FYresp(end,:,4))',s,'b');
    set(h,'EdgeAlpha',.00)
    set(h,'FaceAlpha',.1)
    
    hold off
    set(gca,'XLim',[0 nsteps],'XTick',[0 nsteps/2 nsteps],'FontSize',9)
    title('FFR')
    ss=15;
    % if unique(cell2mat(regexp(varnames(1+zz),'RER')))
    %     ss=16;
    % end
    i=1;
    while i<=ss
        if i+15*(jj-1)<=size(resp,1)
            subplot(4,4,i+1)
            plot(s,[z squeeze(resp(i+zz,:,1))'],'LineWidth',1.5);
            hold on
            h=ciplot(squeeze(resp(i+zz,:,2))',squeeze(resp(i+zz,:,5))',s,'r');
            set(h,'EdgeAlpha',.00)
            set(h,'FaceAlpha',.1)
            h=ciplot(squeeze(resp(i+zz,:,3))',squeeze(resp(i+zz,:,4))',s,'b');
            set(h,'EdgeAlpha',.00)
            set(h,'FaceAlpha',.1)
            set(gca,'XLim',[0 nsteps],'XTick',[0 nsteps/2 nsteps],'FontSize',9)
            title(varnames(i+zz))
            hold off
            i=i+1;
        else
            break
        end
        
    end
    if saveyes==1
        graph2imwrite(fg(jj),['IRF_',num2str(jj)],'png',100,12);
    end
    zz=zz+ss;
end