function [X,y0,k,y1,g1] = litterman(rho,sigma,lambda,varargin)
% litterman  Litterman's prior dummy observations for BVARs.
%
% Syntax
% =======
%
%     [X,Y0,K,Y1,G1] = BVAR.litterman(RHO,SIGMA,LAMBDA,[NY,P])
%     [X,Y0,K,Y1,G1] = BVAR.litterman(RHO,SIGMA,LAMBDA,[NY,P,NG])
%
% Input arguments
% ================
%
% * `RHO` [ numeric ] - White-noise priors (`RHO = 0`) or random-walk
% priors (`RHO = 1`).
%
% * `SIGMA` [ numeric ] - Weight on the dummy observations.
%
% * `LAMBDA` [ numeric ] - Exponential increase in weight depending on the
% lag; `LAMBDA = 0` means all lags are weighted equally.
%
% * `[NY,P,NG]` [ numeric ] - Number of variables, order of the VAR, and
% number of co-integrating vectors in the VAR for which the prior dummies
% will be created.
%
% Output arguments
% =================
%
% * `X` [ numeric ] - Array with prior dummy observations that can be used
% in the `'BVAR='` option of the [`VAR/estimate`](VAR/estimate) function.
%
% * `Y0`, `K` ,`Y1`, `G1` [ numeric ] - These extra output arguments only
% provide more details on the structure of the dummy observations, and have
% no use; `X=[Y0;K;Y1;G1]`.
%
% Description
% ============
%
% Example
% ========
%

% -IRIS Toolbox.
% -Copyright (c) 2007-2012 Jaromir Benes.

ng = 0;
if length(varargin) == 1
   % BVAR.litterman(rho,sigma,lambda,[ny,p]).
   % BVAR.litterman(rho,sigma,lambda,[ny,p,ng]).
   ny = varargin{1}(1);
   p = varargin{1}(2);
   if length(varargin{1}) > 2
      ng = varargin{1}(3);
   end
else
   % Bkw compatibility.
   % BVAR.litterman(rho,sigma,lambda,ny,p).
   % BVAR.litterman(rho,sigma,lambda,ny,p,ng).
   ny = varargin{1};
   p = varargin{2};
   varargin(1:2) = [];
   if length(varargin) > 2
      ng = varargin{3};
   end
end

P = inputParser();
P.addRequired('ny',@isnumericscalar);
P.addRequired('p',@isnumericscalar);
P.addRequired('ng',@isnumericscalar);
P.addRequired('rho',@(x) isnumeric(x) ...
   && (numel(x) == 1 || numel(x) == ny) && all(x >= 0 & x <= 1));
P.addRequired('sigma',@(x) isnumeric(x) ...
   && (numel(x) == 1 || numel(x) == ny) && all(x >= 0));
P.addRequired('lambda',@(x) isnumericscalar(x) && x >= 0);
P.parse(ny,p,ng,rho,sigma,lambda);

%**************************************************************************

rho = rho(:);
sigma = sigma(:);
nsigma = numel(sigma);
if nsigma == 1 && ny > 1
   sigma = sigma(ones(1,ny),1);
elseif nsigma ~= ny
   error('iris:bvar', ...
      ['The number of weights in `sigma` must be scalar ',...
       'or match the number of variables `ny`.']);
end

nd = ny*p;
k = zeros([1,nd]);
g1 = zeros([ng,nd]);
y0 = [diag(sigma.*rho),zeros(ny,nd-ny)];
sigma = sigma(:,ones(1,p));
if lambda > 0
   lags = (1 : p).^lambda;
   lags = lags(ones(1,ny),:);
   sigma = sigma .* lags;
end
y1 = diag(sigma(:));  
X = [y0;k;y1;g1];

end