% frbsffavar_China.m
%
% Forecast a wide variety of macroeconomic data series using a dynamic factor
%  model/factor-augmented VAR (FAVAR).  This is version 1.0.
%
% written by Eric T. Swanson, 2006.
% modified by Kuni Natsuki

clear all ;


%%%%%%%%%%%%Input file name%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Name of input data file; default is 'Input.xls'
infilename = 'Input.xls';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Here are the key parameters of the FAVAR/Dynamic Factor Model (edit as desired):
%

infiledir = '' ;  % name of input data file directory; default is 'data'
outfiledir = '' ; % name of output results file directory; default is 'results/(today's date)'
Neafactors = 2 ;  % number of Economic Activity factors to extract
Nprfactors = 0 ;  % number of Inflation factors to extract
Nintfactors = 0 ; % number of Interest Rate Factors to extract
Ntrafactors = 0 ; % number of International Trade Factors to extract
varlags = 12 ; % number of monthly lags to use in the FAVAR for forecasting
fcasthorizon = 36 ; % number of months ahead to forecast
factorstartdate = [] ; % [year month] to begin the factor extraction.
               %  Note: means and HP trends (if any) will still be computed from earliest
               %  possible month, regardless of factorstartdate.
varstartdate = [2000 5] ; % [year month] to begin the estimation of the VAR (FAVAR).
               %  Note: the factor extraction is peformed on the data back to
               %  factorstartdate.  varstartdate must be >= factorstartdate.
enddate = [] ; % [year month] to end the factor extraction and the FAVAR estimation
               %  (this is useful for computing historical out-of-sample forecasts).
               %  [] uses data up through the current month.  Note: means and HP trends
               %  (if any) will ignore all data after enddate (for out-of-sample forecasting).
diffspecification = 'no' ; % if 'yes', then all HP detrending will be abandoned and a
               %  3-month log-first difference specification will be used for those
               %  variables instead.

%
% Don't edit the code from here on unless you really know what you're doing!!
%

headrows = 3 ; % number of non-numeric rows at top of haverseries Excel file
burn = 3 ; % number of observations at beginning of data to "burn".  Burn must be
            %  >= number of months of differencing; i.e., if there are any annual
            %  difference transformations being done to the data, then burn must
            %  be >=12.  Note: if factorstartdate has been set, then burn is automatically
            %  updated below to factorstartdate - first date of data.

fprintf('\n\nRunning FRBSF Dynamic Factor Model/FAVAR version 1.0:\n\n') ;

% Read the data in from Excel:
func_dir=[cd '\library'];
addpath(func_dir);

if strcmp(infiledir,''); infiledir = 'data'; end ;
if strcmp(infilename,''); infilename = 'Input.xls'; end ;
infilestr = strcat(infiledir,'/',infilename) ;
fprintf('Reading data from Excel spreadsheet: %s ...\n',infilestr) ;
[data,txt,raw] = xlsread(infilestr,'Monthly') ;
dates = str2num(cell2mat(txt(headrows+1:size(txt,1),1))) ;
dates = [floor(dates/100) mod(dates,100)] ;

% Create directory for favar forecast, data archive, and log files:
outfiledir = strcat('results/',datestr(now,'yyyy.mm.dd')); 

if exist(outfiledir, 'dir'); % if, say, results\20130702 exists, then create results\20130702.1
        run_number = 1;
        while exist([outfiledir '.' num2str(run_number) ], 'dir')
            run_number = 1 + run_number;
        end % while
        this_dir = [outfiledir '.' num2str(run_number) ];
else
        this_dir = outfiledir;
end
 outfiledir = this_dir;



fprintf('Ouput will be written to directory: %s\n',outfiledir) ;
warning off MATLAB:MKDIR:DirectoryExists ;
mkdir(outfiledir); mkdir(strcat(outfiledir,'/figures')); mkdir(strcat(outfiledir,'/figuresdiff')) ;
outfilestr = strcat(outfiledir,'/favarfcast.', datestr(now,'yyyy.mm.dd')) ;
if strcmp(diffspecification,'yes'); outfilestr = strcat(outfilestr,'diff.log') ;
else outfilestr = strcat(outfilestr,'.log') ;
end ;
fprintf('Terminal output being logged to file: %s\n',outfilestr) ;
warning off MATLAB:DELETE:FileNotFound ;
delete(outfilestr) ; % if a diary file of this name exists already, delete it
diary(outfilestr) ;

% Record parameter values to log file:
fprintf('Key Parameters of FAVAR are as Follows\n') ;
fprintf('Number of Economic Activity Factors   :  %i\n', Neafactors) ;
fprintf('Number of Inflation Factors           :  %i\n', Nprfactors) ;
fprintf('Number of Interest Rate Factors       :  %i\n', Nintfactors) ;
fprintf('Number of International Trade Factors :  %i\n', Ntrafactors) ;
fprintf('Factor Extraction Start Date          :  ') ;
if(isempty(factorstartdate)); fprintf('[]\n'); else; fprintf('[%i %i]\n', factorstartdate); end ;
fprintf('VAR (FAVAR) Estimation Start Date     :  ') ;
if(isempty(varstartdate)); fprintf('[]\n'); else; fprintf('[%i %i]\n', varstartdate); end ;
fprintf('Factor Extraction and VAR End Date    :  ') ;
if(isempty(enddate)); fprintf('[]\n'); else; fprintf('[%i %i]\n', enddate); end ;
fprintf('Number of (monthly) Lags in VAR       :  %i\n', varlags) ;
fprintf('Forecast Horizon (in months)          :  %i\n', fcasthorizon) ;
fprintf('No-Gaps First-Difference Specification:  %s\n', diffspecification) ;

% Impose factorstartdate and enddate
% If factorstartdate is set, then that tells us how many observations to burn (after detrending):
if (~isempty(factorstartdate)) ;
  burn = 12*(factorstartdate(1)-dates(1,1)) + factorstartdate(2)-dates(1,2) ;
end ;
if (~isempty(enddate)) ; % chop off data at enddate:
  i = find(dates(:,1)==enddate(1) & dates(:,2)==enddate(2)) ;
  data = data(1:i,:) ;
  txt = txt(1:headrows+i,:) ;
  raw = raw(1:headrows+i,:) ;
  dates = dates(1:i,:) ;
end ;

Nseries = size(data,2) ;
T = size(data,1) - burn ;
for i = 1:Nseries ;
  j = find(~isnan(data(burn+1:burn+T,i))) ;
  firstobs(i) = j(1) ; % first observation for each column of data
  lastobs(i) = j(end) ; % last observation for each column of data
end ;

%
% Step 1: Transform the data by detrending or first-differencing each group of
%  series.  Then normalize each series to have mean zero and unit standard deviation.
%  Hold on to the trends, means, and standard deviations so that we can transform
%  the data back to its original format at the very end.
%

fprintf('Standardizing data for factor analysis...\n') ;

% Divide the data series into various transformations to be performed:
statindexes = find(strcmp(raw(2,3:end),'NONE')) ;
hplindexes = find(strcmp(raw(2,3:end),'HPLOOSE')) ;
hptindexes = find(strcmp(raw(2,3:end),'HPTIGHT')) ;
diffindexes = find(strcmp(raw(2,3:end),'DIFF')) ;
normindexes = find(strncmp(raw(2,3:end),'NORM',4)) ;

if strcmp(diffspecification,'yes') ; % if difference specification is requested
  diffindexes = [diffindexes hplindexes hptindexes] ;
  difflevindexes = normindexes ; % use level (not log) diff on normalized series
  hplindexes = [] ;
  hptindexes = [] ;
else ;
  hptindexes = [hptindexes,normindexes] ; % use HPTIGHT on normalized series
  difflevindexes = [] ;
end ;

% normalize series that require normalization (inventories, trade):
for i = 1:length(normindexes) ;
  normby = char(raw(2,2+normindexes(i))) ; % Haver mnemonic for series to normalize by
  normbyindexes(i) = find(strcmp(raw(3,3:end),normby(6:end-1))) ;
end ;
%data(:,normindexes) = data(:,normindexes) ./ data(:,normbyindexes);

% no transformation: just demean
statmean = robustmean(data(:,statindexes)) ; % mean
datastat = data(burn+1:burn+T,statindexes) - repmat(statmean,T,1) ; % demean

% HP loose: this is for data that show a strong trend, thus we log the data first,
%  then use HP parameter of 5e5; the 3e7 is for longer-term forecasting
trendhpl = hpfit(log(data(:,hplindexes)),5e5);  trendhpl = trendhpl(burn+1:burn+T,:) ;
trendhpllt = hpfit(log(data(:,hplindexes)),3e7); trendhpllt = trendhpllt(burn+1:burn+T,:) ;
datahpl = log(data(burn+1:burn+T,hplindexes)) - trendhpl ; % detrend

% HP tight: this is for data that do not have a strong trend;
%  use HP parameter of 5e7; the 3e8 is for longer-term forecasting
trendhpt = hpfit(data(:,hptindexes),5e7); trendhpt = trendhpt(burn+1:burn+T,:) ;
trendhptlt = hpfit(data(:,hptindexes),3e8); trendhptlt = trendhptlt(burn+1:burn+T,:) ;
datahpt = data(burn+1:burn+T,hptindexes) - trendhpt ; % detrend

% difference n: take log-difference of the data from n months previous:
datadiff = log(data(burn+1:burn+T,diffindexes)) - log(data(burn-2:end-3,diffindexes)) ;
diffmean = robustmean(datadiff) ;
datadiff = datadiff - repmat(diffmean,T,1) ;

% level difference: take level (not log) diff of the data from n months previous:
datadifflev = data(burn+1:burn+T,difflevindexes) - data(burn-2:end-3,difflevindexes) ;
difflevmean = robustmean(datadifflev) ;
datadifflev = datadifflev - repmat(difflevmean,T,1) ;

% reassemble data set and standardize to unit standard deviation:
datatr(:,[statindexes,hplindexes,hptindexes,diffindexes,difflevindexes]) = ...
					[datastat,datahpl,datahpt,datadiff,datadifflev] ;
datastd = robuststd(datatr) ;
datatr = datatr ./ repmat(datastd,T,1) ;

% finally, divide the data series into groups (or categories) of series for factor analysis:
eaindexes = find(strcmp(raw(1,3:end),'EA')) ;
prindexes = find(strcmp(raw(1,3:end),'PR')) ;
intindexes = find(strcmp(raw(1,3:end),'INT')) ;
traindexes = find(strcmp(raw(1,3:end),'TRA')) ;

econactivity = datatr(:,eaindexes) ;
prices = datatr(:,prindexes) ;
intrates = datatr(:,intindexes) ;
trade = datatr(:,traindexes) ;


%
% Step 2: In general, there will be some series for which we do not yet have
%  observations at the end of the sample.  We fill in these missing values with
%  optimal forecasts based on the model.
% In periods 1 through T-2, missing data is sparse, so we fill in those missing
%  data using static principal components with the function fillinmissdata.m
% In periods T-1 and T, it is often the observed data rather than the missing data that
%  is sparse.  Thus, we will use the dynamic structure of the model to help estimate
%  the factors in those periods.  This is done with the function dynfactorfcast.m
%

% Fill in missing data up through min(lastobs) using static principal components:
T1 = min(lastobs) - 1 ;
fprintf('Filling in missing EA data:\n') ;
 econactivity(1:T1,:) = favarfillindata(econactivity(1:T1,:), Neafactors) ;
fprintf('Filling in missing PR data:\n') ;
 prices(1:T1,:) = favarfillindata(prices(1:T1,:), Nprfactors) ;
fprintf('Filling in missing INT data:\n') ;
 intrates(1:T1,:) = favarfillindata(intrates(1:T1,:), Nintfactors) ;
fprintf('Filling in missing TRA data:\n') ;
 trade(1:T1,:) = favarfillindata(trade(1:T1,:), Ntrafactors) ;
 
% Now fill in data for the most recent periods (the more interesting case):
fprintf('\nFilling in data in most recent periods using MLE/Kalman filter if possible:\n') ;
idiosyncerrwhite = repmat(NaN,T-1,Nseries) ; % initialize idiosyncerrwhite
for T1 = min(lastobs):T ;
  % To fill in missing data at time T1, we extract factors up through time T1-1:
  [eafactors,ealoadings,dum,eaerrors] = extract(econactivity(1:T1-1,:), Neafactors) ;
  [prfactors,prloadings,dum,prerrors] = extract(prices(1:T1-1,:), Nprfactors) ;
  [intfactors,intloadings,dum,interrors] = extract(intrates(1:T1-1,:), Nintfactors) ;
  [trafactors,traloadings,dum,traerrors] = extract(trade(1:T1-1,:), Ntrafactors) ;
  idiosyncerrors(:,[eaindexes,prindexes,intindexes,traindexes]) = ...
                                         [eaerrors,prerrors,interrors,traerrors] ;

  % then estimate the dynamic factor model up through T1-1 using a VAR (FAVAR):
  vardata = [eafactors,prfactors,intfactors,trafactors] ;
  [B,dum,dum2,dum3,omegaeps] = VAR(vardata,varlags,0) ;

  % then estimate the serial correlation of the idiosyncratic errors:
  for i = find(lastobs>=T1) ;
    [rho(i),dum,dum2,idiosyncerrwhite(firstobs(i)+1:T1-1,i)] = ...
        ols(idiosyncerrors(firstobs(i)+1:T1-1,i),idiosyncerrors(firstobs(i):T1-2,i),0) ;
  end ; % serial correlation is saved in rho
  idiosynccov = robustcov(idiosyncerrwhite) ; % covariance matrix of whitened idiosync errors

  % compute the forecast for period T1 implied by the VAR (FAVAR):
  if (~isempty(varstartdate)); j = 12*(varstartdate(1)-dates(1,1)) + varstartdate(2)-dates(1,2) - burn ;
  else; j = 1 ;
  end ;
  varfcast = VARfcast(vardata(j:end,:),varlags,1) ;
  eavarind = [1:Neafactors] ;
  prvarind = Neafactors + [1:Nprfactors] ;
  intvarind = Neafactors + Nprfactors + [1:Nintfactors] ;
  travarind = Neafactors + Nprfactors + Nintfactors + [1:Ntrafactors] ;

  % combine the VAR forecast with the observed time-T1 data to estimate the
  %  factors at time T1 (using MLE/Bayesian updating/Kalman filter if possible):
  fprintf('Filling in missing data at time T'); if (T1<T); fprintf('%i',T1-T); end; fprintf('...\n');
  fprintf(' estimating EA factors at time T'); if (T1<T); fprintf('%i',T1-T); end; fprintf(':\n');
  eaobs = ~isnan(econactivity(T1,:)) ; % locations of the non-missing data
  eaobs2 = eaindexes(eaobs) ;
  eafactors(T1,:) = favardynfactfcast(varfcast(:,eavarind),omegaeps(eavarind,eavarind),...
             econactivity(T1,eaobs),ealoadings(:,eaobs),...
             rho(eaobs2).*idiosyncerrors(T1-1,eaobs2),idiosynccov(eaobs2,eaobs2)) ;
  fprintf(' estimating PR factors at time T'); if (T1<T); fprintf('%i',T1-T); end; fprintf(':\n');
  probs = ~isnan(prices(T1,:)) ;
  probs2 = prindexes(probs) ;
  prfactors(T1,:) = favardynfactfcast(varfcast(:,prvarind),omegaeps(prvarind,prvarind),...
             prices(T1,probs),prloadings(:,probs),...
             rho(probs2).*idiosyncerrors(T1-1,probs2),idiosynccov(probs2,probs2)) ;
  fprintf(' estimating INT factors at time T'); if (T1<T); fprintf('%i',T1-T); end; fprintf(':\n');
  intobs = ~isnan(intrates(T1,:)) ;
  intobs2 = intindexes(intobs) ;
  intfactors(T1,:) = favardynfactfcast(varfcast(:,intvarind),omegaeps(intvarind,intvarind),...
             intrates(T1,intobs),intloadings(:,intobs),...
             rho(intobs2).*idiosyncerrors(T1-1,intobs2),idiosynccov(intobs2,intobs2)) ;
  fprintf(' estimating TRA factors at time T'); if (T1<T); fprintf('%i',T1-T); end; fprintf(':\n');
  traobs = ~isnan(trade(T1,:)) ;
  traobs2 = traindexes(traobs) ;
  trafactors(T1,:) = favardynfactfcast(varfcast(:,travarind),omegaeps(travarind,travarind),...
             trade(T1,traobs),traloadings(:,traobs),...
             rho(traobs2).*idiosyncerrors(T1-1,traobs2),idiosynccov(traobs2,traobs2)) ;

  % finally, infer missing data and missing idiosyncratic errors at time T1:
  idiosyncerrors(T1,:) = idiosyncerrors(T1-1,:) .* rho ;
  econactivity(T1,~eaobs) = eafactors(T1,:) * ealoadings(:,~eaobs) + idiosyncerrors(T1,eaindexes(~eaobs)) ;
  prices(T1,~probs) = prfactors(T1,:) * prloadings(:,~probs) + idiosyncerrors(T1,prindexes(~probs)) ;
  intrates(T1,~intobs) = intfactors(T1,:) * intloadings(:,~intobs) + idiosyncerrors(T1,intindexes(~intobs)) ;
  trade(T1,~traobs) = trafactors(T1,:) * traloadings(:,~traobs) + idiosyncerrors(T1,traindexes(~traobs)) ;
end ; % end for T1 = min(lastobs):T

idiosyncerrors(T,eaobs2) = econactivity(T,eaobs) - eafactors(T,:) * ealoadings(:,eaobs) ;
idiosyncerrors(T,probs2) = prices(T,probs) - prfactors(T,:) * prloadings(:,probs) ;
idiosyncerrors(T,intobs2) = intrates(T,intobs) - intfactors(T,:) * intloadings(:,intobs) ;
idiosyncerrors(T,traobs2) = trade(T,traobs) - trafactors(T,:) * traloadings(:,traobs) ;

%
% Step 3: Forecast the factors using the dynamic factor model.
%

fprintf('\nUsing VAR (FAVAR) to forecast factors out to %i months...\n',fcasthorizon) ;
vardata = [eafactors,prfactors,intfactors,trafactors] ;
varfcast = VARfcast(vardata,varlags,fcasthorizon,B) ;
eafactors(T+1:T+fcasthorizon,:) = varfcast(:,eavarind) ;
prfactors(T+1:T+fcasthorizon,:) = varfcast(:,prvarind) ;
intfactors(T+1:T+fcasthorizon,:) = varfcast(:,intvarind) ;
trafactors(T+1:T+fcasthorizon,:) = varfcast(:,travarind) ;

%
% Step 4: Forecast all of the observable variables.
%  In order to forecast all of the observable variables, we need to do two things:
%  First, we need to forecast the idiosyncratic error terms for each observable series.
%   If the last observation of an observed variable is very different from what
%   the factor model would predict, then we want to forecast that residual to
%   persist for some time rather than to die out immediately in the next period.
%  Second, we need to transform observed variables back to their non-zero-mean,
%   non-unit-std-deviation versions.  The hardest part of this is projecting the
%   HP trends forward to the end of the forecast horizon.
%

% forecast idiosyncratic errors using serial correlation in the idiosyncratic error terms:
fprintf('Forecasting idiosyncratic errors (using serial correlation)...\n') ;
idiosyncerrors(T+1:T+fcasthorizon,:) = repmat(idiosyncerrors(T,:),fcasthorizon,1) ...
               .* repmat(rho,fcasthorizon,1).^repmat([1:fcasthorizon]',1,Nseries) ;

% forecast all observable variables:
fprintf('Forecasting all observable variables...\n') ;
econactivity(T+1:T+fcasthorizon,:) = eafactors(T+1:T+fcasthorizon,:) * ealoadings + ...
                                          idiosyncerrors(T+1:T+fcasthorizon,eaindexes) ;
prices(T+1:T+fcasthorizon,:) = prfactors(T+1:T+fcasthorizon,:) * prloadings + ...
                                          idiosyncerrors(T+1:T+fcasthorizon,prindexes) ;
intrates(T+1:T+fcasthorizon,:) = intfactors(T+1:T+fcasthorizon,:) * intloadings + ...
                                          idiosyncerrors(T+1:T+fcasthorizon,intindexes) ;
trade(T+1:T+fcasthorizon,:) = trafactors(T+1:T+fcasthorizon,:) * traloadings + ...
                                          idiosyncerrors(T+1:T+fcasthorizon,traindexes) ;

% Transform observed variables back to non-zero-mean, non-unit-std-dev versions:
datatr(1:T+fcasthorizon,[eaindexes,prindexes,intindexes,traindexes]) = ...
                                          [econactivity,prices,intrates,trade] ;
datatr = datatr .* repmat(datastd,T+fcasthorizon,1) ;

% no transformation: just restore the mean
datastat = datatr(:,statindexes) + repmat(statmean,T+fcasthorizon,1) ;

% HP loose: forecast and restore the HP trend
trendhpl = trendfcast(trendhpl,fcasthorizon,trendhpllt) ;
datahpl = exp(datatr(:,hplindexes) + trendhpl) ;

% HP tight: forecast and restore the HP trend
trendhpt = trendfcast(trendhpt,fcasthorizon,trendhptlt) ;
datahpt = datatr(:,hptindexes) + trendhpt ;

% diff: restore mean and transform data back to levels
datadiff = datatr(:,diffindexes) + repmat(diffmean,T+fcasthorizon,1) ;
for i = 1:length(diffindexes) ;
  T1 = find(isfinite(data(burn+1:end,diffindexes(i)))) ;
  datalevel(T1(1):T1(1)+2,i) = data(T1(1)+burn:T1(1)+burn+2,diffindexes(i)) ;
  for T2 = T1(1)+3:T+fcasthorizon ; % forecast level forward
    datalevel(T2,i) = datalevel(T2-3,i) * exp(datadiff(T2,i)) ;
  end ;
  for T2 = T1(1)-1:-1:1 ; % forecast level backward
    datalevel(T2,i) = datalevel(T2+3,i) / exp(datadiff(T2+3,i)) ;
  end ;
end ;
% datadiff = datalevel; clear datalevel ;
% 
% % diff: restore mean and transform data back to levels
% datadifflev = datatr(:,difflevindexes) + repmat(difflevmean,T+fcasthorizon,1) ;
% for i = 1:length(difflevindexes) ;
%   T1 = find(isfinite(data(burn+1:end,difflevindexes(i)))) ;
%   datalevel(T1(1):T1(1)+2,i) = data(T1(1)+burn:T1(1)+burn+2,difflevindexes(i)) ;
%   for T2 = T1(1)+3:T+fcasthorizon ; % forecast level forward
%     datalevel(T2,i) = datalevel(T2-3,i) + datadifflev(T2,i) ;
%   end ;
%   for T2 = T1(1)-1:-1:1 ; % forecast level backward
%     datalevel(T2,i) = datalevel(T2+3,i) - datadifflev(T2+3,i) ;
%   end ;
% end ;
% if (~isempty(difflevindexes)); datadifflev = datalevel; clear datalevel ;
% else; datadifflev = repmat(NaN,T+fcasthorizon,0) ;
% end ;
    
% recombine data into matrix form:
datatr(:,[statindexes,hplindexes,hptindexes,diffindexes,difflevindexes]) = ...
            					[datastat,datahpl,datahpt,datadiff,datadifflev] ;

% un-normalize series that required normalization (inventories, trade):
%datatr(:,normindexes) = datatr(:,normindexes) .* datatr(:,normbyindexes) ;

eafinal = datatr(:,eaindexes) ; % useful to save these series for graphing
prfinal = datatr(:,prindexes) ;
intfinal = datatr(:,intindexes) ;
trafinal = datatr(:,traindexes) ;

%
% Step 5: Record the input data and forecast data to an Excel workbook with a date
%  stamp for posterity and for possible graphing, analysis, etc.
% The first worksheet in the workbook reports the input data (including missing data)
% The second worksheet in the workbook reports the forecasts for the observabe data 
%  (including history), with all missing data at end of sample filled in
% The third worksheet reports the factor estimates and factor forecasts
%

fprintf('Writing data to Excel workbook...\n') ;
outfilestr(end-2:end) = 'xls' ; % replace .log extension with .xls
copyfile('data/template_temp.xls',outfilestr) ;

% fix glitch in Excel/Matlab #N/A interface by replacing error strings with #N/A:
raw(find(strcmp(raw,'ActiveX VT_ERROR: '))) = cellstr('#N/A') ;

% read in dates and extend them forward to end of forecast horizon:
dates(burn+T+1:burn+T+fcasthorizon,2) = mod(dates(burn+T,2) +[1:fcasthorizon]'-1,12) +1 ;
dates(burn+T+1:burn+T+fcasthorizon,1) = dates(burn+T,1) + floor((dates(burn+T,2)+[1:fcasthorizon]'-1)/12) ;
datesout1 = [100*dates(:,1) + dates(:,2)] ;
datesout2 = reshape(sprintf('%02i/1/%i', (fliplr(dates))'),9,burn+T+fcasthorizon)' ;

% output observable data history and forecasts to second worksheet of workbook:
xlswrite(outfilestr,raw,'Data and Forecasts') ;
xlswrite(outfilestr,num2cell(datatr),'Data and Forecasts',sprintf('c%i',headrows+burn+1)) ;
xlswrite(outfilestr,num2cell(datesout1),'Data and Forecasts',sprintf('a%i',headrows+1)) ;
xlswrite(outfilestr,cellstr(datesout2),'Data and Forecasts',sprintf('b%i',headrows+1)) ;

% output factor estimates and forecasts to third worksheet of workbook:
for i=1:Neafactors; titlestr{i} = sprintf('EA%i',i); end ;
for i=1:Nprfactors; titlestr{Neafactors+i} = sprintf('PR%i',i); end ;
for i=1:Nintfactors; titlestr{Neafactors+Nprfactors+i} = sprintf('INT%i',i); end ;
for i=1:Ntrafactors; titlestr{Neafactors+Nprfactors+Nintfactors+i} = sprintf('TRA%i',i); end ;
titlestr(2:burn+1,:) = cellstr('#N/A') ;
xlswrite(outfilestr,titlestr,'Factors','b1') ;
xlswrite(outfilestr,cellstr(datesout2),'Factors','a2') ;
xlswrite(outfilestr,num2cell([vardata;varfcast]),'Factors',sprintf('b%i',2+burn)) ;

% output quarterly forecasts to fourth worksheet of workbook:
fprintf('Forecasting GDP...\n') ;
%favarfcastgdp_121227 ; % call favarfcastgdp script
fprintf('Creating FRBSF Quarterly Forecast Sheet...\n') ;

xlswrite(outfilestr,cellstr(datesout2),'GDP','a2') ;
xlswrite(outfilestr,cellstr('Chinese GDP'),'GDP','b1') ;
%xlswrite(outfilestr,num2cell(gdpm),'GDP',sprintf('b%i',2+burn)) ;

%Print out the loadings
xlswrite(outfilestr,raw(3,(3:eaindexes(end)+2)),'Loadings','b1') ;
for j=1:Neafactors
xlswrite(outfilestr,cellstr(sprintf('EA%i',j)),'Loadings',sprintf('a%i',j+1)) ;
xlswrite(outfilestr,num2cell(ealoadings),'Loadings',sprintf('b%i',j+1)) ;

    %What are the factors' variance shares?
    mean(ealoadings.^2,2)
end

if Nprfactors > 0
xlswrite(outfilestr,raw(3,(eaindexes(end)+3):(prindexes(end)+2)),'Loadings',sprintf('b%i',j+1+1+1)) ;    
for k=1:Nprfactors
    xlswrite(outfilestr,cellstr(sprintf('PR%i',k)),'Loadings',sprintf('a%i',j+1+1+1+k)) ;
    xlswrite(outfilestr,num2cell(prloadings),'Loadings',sprintf('b%i',j+1+1+1+k)) ;
end
end

if Nintfactors >0
xlswrite(outfilestr,raw(3,(prindexes(end)+3):(intindexes(end)+2)),'Loadings',sprintf('b%i',j+1+1+1+k+1+1)) ; 
for l=1:Nintfactors
    xlswrite(outfilestr,cellstr(sprintf('INT%i',l)),'Loadings',sprintf('a%i',j+1+1+1+k+1+1+l)) ;
    xlswrite(outfilestr,num2cell(intloadings),'Loadings',sprintf('b%i',j+1+1+1+k+1+1+l)) ;
end    
    
end



%Not used for Chinese FAVAR
%frbsffcastsheet ; % call frbsffcastsheet script

% archive program files and workspace to directory as well
fprintf('Archiving data to %s\n',outfiledir) ;
copyfile(infilestr,outfiledir) ;

fprintf('Archiving FAVAR code to %s\n',outfiledir) ;
copyfile(func_dir,[outfiledir '\library']) ;

date = dates(burn+1:end,1) + dates(burn+1:end,2)/12 - 1/24 ; % date vector for plotting

fprintf('Archiving Matlab workspace to %s\n',outfiledir) ;
save(outfilestr(1:end-4)) ;

fprintf('Done.\n\n') ;
diary off ; % turn off log file
warning on MATLAB:MKDIR:DirectoryExists ;
warning on MATLAB:DELETE:FileNotFound ; % reset warning status





