function fcast = favardynfactfcast(varfcast,varfcastcov,observs,loadings,idiosyncerrs,idiosynccov) ;
% function fcast = favardynfactfcast(varfcast,omegaeps,observs,loadings,idiosyncerrs,idiosynccov) ;
%
% In a dynamic factor model with data being released sequentially in real time,
%  there will typically be missing data at the end of the sample.  This routine
%  fills in those missing data with optimal (or near optimal) forecasts.
% The function considers a static forecast (based only on data at time T) and a
%  dynamic forecast (based on data through time T-1, projected into period T
%  using a VAR).  Ideally, we would combine both forecasts optimally to forecast
%  the data in period T.  However, doing WLS/GLS on the time-T observable data
%  will often be infeasible in practice due to the enormous number of observable
%  series that are typically present in a dynamic factor model.  At some point
%  (for example, when n > T), the estimated variance matrix of the observables will
%  surely be singular, and for much smaller values of n the variance matrix will
%  be very poorly estimated and non-robust despite being technically nonsingular.
% Thus, the routine will only do WLS/GLS on the time-T observable data only for
%  relatively small numbers of non-missing time-T observations.  Beyond that point,
%  the routine will switch over to pure OLS on the time-T observations alone to
%  estimate the time-T factor value.  This estimate is consistent in n and is
%  going to be much more numerically robust for moderate to large n.
% The function input varfcast is a dynamic VAR/FAVAR forecast based on data up
%  through time T-1, treating the factors up through time T-1 as observed.
% omegaeps is the variance matrix of the FAVAR shocks/residuals.  It thus gives
%  us the uncertainty surrounding the dynamic forecast.
% observables is the vector of time-T observations.
% loadings is the matrix of factor loadings of the observable variables on the
%  factors.
% errorvar is the variance matrix of the idiosyncratic component of the observable
%  variables.  It thus provides a measure of uncertainty surrounding the static
%  forecast (based on the time-T observed data).


k = size(varfcast,2) ; % length of (row) vector to forecast
n = size(observs,2) ; % number of period-T observations

if (n==0) ; % if no time-T observables, estimate factor dynamically (using VAR):
  fprintf('  no current period observations, inferring factors from last period observations...\n') ;
  fcast = varfcast ;
elseif (cond(idiosynccov)>10000) ; % if variance matrix of observables is poorly
                              % conditioned, then get factor by OLS, which is
                              % unbiased and consistent in the number of observables.
                              % Can also try elseif (cond(errorvar)>10000 | n>2*k) ;
  fprintf('  condition number %6.0f, inferring factors using OLS...\n',cond(idiosynccov)) ;
  fcast = ols(observs',loadings',0)' ;
else ; % if only a few observables, then estimate factor jointly by
       %  MLE/Bayesian updating/Kalman Filter:
  fprintf('  condition number %6.1f, inferring factors using MLE/Kalman filter...\n',cond(idiosynccov)) ;
  fcast = varfcast + (observs -varfcast*loadings -idiosyncerrs) / ...
                        (loadings'*varfcastcov*loadings + idiosynccov) *loadings'*varfcastcov ;
% The above formula is essentially Bayes' rule. Here's an alternative formula based on maximum
%  likelihood which is mathematically identical but probably less robust numerically, due to
%  extra matrix inversions (would need to incorporate idiosyncerrs into this formula):
%  fcast = (varfcast /varfcastcov + observs /idiosynccov * loadings') ...
%                         / (inv(varfcastcov) + loadings /idiosynccov * loadings') ;
end ;



