%
% Compute GDP, output gap, and potential GDP forecast:
%


% load GDP data
[dataq,txtq,rawq] = xlsread(infilestr,'Quarterly') ;
datesq = str2num(cell2mat(txtq(headrows+1:size(txtq,1),1))) ;
datesq = [floor(datesq/10) mod(datesq,10)] ;
%burnq = find(datesq(:,1)==dates(burn+1,1) & datesq(:,2)==floor(dates(burn,2)/3)) ;
burnq=1;
Tq = find(~isnan(dataq)); Tq = Tq(end) - burnq ; % last GDP observation
firstm = find(mod(dates(burn+1:burn+3,2),3)==1) ; % first month we'll impute monthly GDP
lastm = find(dates(:,1)==datesq(Tq+burnq,1) & dates(:,2)==datesq(Tq+burnq,2)*3) - burn ; % last month we'll impute monthly GDP

% detrend (or difference) GDP data:
if strcmp(diffspecification,'yes') ;
  Dlgdp = log(dataq(burnq+1:burnq+Tq)) - log(dataq(burnq:burnq+Tq-1)) ; % 1-quarter difference
  Dlgdpmean = robustmean(Dlgdp) ;
  gdpstd = Dlgdp - repmat(Dlgdpmean,Tq,1) ;
else ;
  lgdp = log(dataq) ;
  trendgdp = hpfit(lgdp); trendgdp = trendgdp(burnq+1:burnq+Tq) ;
  trendgdplt = hpfit(lgdp,1e6); trendgdplt = trendgdplt(burnq+1:burnq+Tq) ;
  gdpstd = lgdp(burnq+1:burnq+Tq) - trendgdp ;
end ;

% estimate GDP loadings on EA factors (which implies a monthly GDP forecast):
eafactorsq = squeeze(mean(reshape(eafactors(firstm:firstm+Tq*3-1,:),3,Tq,Neafactors))) ;

%Convert growth to levels
eafactorsq_length = length(eafactorsq);
eafactorsq_level = zeros(eafactorsq_length,1);
eafactorsq_level(1)=(1+eafactorsq(1)/100)^(1/4);
for i=2:eafactorsq_length
    eafactorsq_level(i)=eafactorsq_level(i-1)*(1+eafactorsq(i)/100)^(1/4);
end

%[gdploadings,dum,dum2,resids] = ols(gdpstd,eafactorsq,0) ;
[gdploadings,dum,dum2,resids] = ols(gdpstd,eafactorsq_level,0) ;

% forecast GDP residuals using serial correlation:
rhogdp = ols(resids(2:Tq),resids(1:Tq-1),0) ;

fcasthorizon=fcasthorizon+36;

resids(Tq+1:Tq+3+ceil(fcasthorizon/3)) = resids(Tq) * rhogdp.^[1:3+ceil(fcasthorizon/3)] ;

fcasthorizon=fcasthorizon-36;
residsm = zeros(T+fcasthorizon,1) ;

for i = firstm:T+fcasthorizon ;
  residsm(i) = resids(ceil((i-firstm+1)/3)) ;
end ;

% compute monthly standardized GDP using gdploadings:
gdpstdm = eafactors * gdploadings + residsm ;

% transform GDP data back to monthly levels and monthly forecasts:
if strcmp(diffspecification,'yes') ; % restore mean and transform data back to levels
  gdpstdm = (gdpstdm + repmat(Dlgdpmean,T+fcasthorizon,1)) /3 ;
  lgdpm = repmat(NaN,T+fcasthorizon,1) ; % initialize monthly gdp variable
  lgdpm(firstm+Tq*3-2) = log(dataq(burnq+Tq)) - 1/3 *gdpstdm(firstm+Tq*3-2) ...
           + 1/3 *gdpstdm(firstm+Tq*3-1) ; % benchmark last quarter to quarterly data
  for i = firstm+Tq*3-1:T+fcasthorizon ; % forecast forward to end of fcasthorizon
    lgdpm(i) = lgdpm(i-1) + gdpstdm(i) ;
  end ;
  for i = firstm+Tq*3-3:-1:firstm ; % forecast backward to beginning of data
    lgdpm(i) = lgdpm(i+1) - gdpstdm(i+1) ;
  end ;
else ; % interpolate GDP trend to monthly and project trend forward
  gdpgapm = gdpstdm ;
  trendgdpm = repmat(NaN,lastm,1) ;
  trendgdpm(firstm:lastm,1) = trendinterp(trendgdp) ;
  trendgdpltm(firstm:lastm,1) = trendinterp(trendgdplt) ;
  trendgdpm = trendfcast(trendgdpm,T+fcasthorizon-lastm,trendgdpltm) ;
  lgdpm = gdpgapm + trendgdpm ;
end ;

gdpm = exp(lgdpm) ;


