function [By,Bx,u,xxi]=rfvar3(ydata,lags,xdata,breaks,lambda,mu)
%function [By,Bx,u,xxi]=rfvar3(ydata,lags,xdata,breaks,lambda,mu)
% This algorithm goes for accuracy without worrying about memory requirements.
% ydata:   dependent variable data matrix
% xdata:   exogenous variable data matrix
% lags:    number of lags
% breaks:  rows in ydata and xdata after which there is a break.  This allows for
%          discontinuities in the data (e.g. war years) and for the possibility of
%          adding dummy observations to implement a prior.  This must be a column vector.
%          Note that a single dummy observation becomes lags+1 rows of the data matrix,
%          with a break separating it from the rest of the data.  The function treats the 
%          first lags observations at the top and after each "break" in ydata and xdata as
%          initial conditions. 
% lambda:  weight on "co-persistence" prior dummy observations.  This expresses
%          belief that when data on *all* y's are stable at their initial levels, they will
%          tend to persist at that level.  lambda=5 is a reasonable first try.
% mu:      weight on "own persistence" prior dummy observation.  Expresses belief
%          that when y_i has been stable at its initial level, it will tend to persist
%          at that level, regardless of the values of other variables.  There is
%          one of these for each variable.  A reasonable first guess is mu=2.
%      Both mu and lambda dummies include coefficients on x.  If you want to exclude
%      them you will have to construct your own dummies, using breaks.
%      The program assumes that the first lags rows of ydata and xdata are real data, not dummies.
%      Dummy observations should go at the end, if any.  If pre-sample x's are not available,
%      repeating the initial xdata(lags+1,:) row or copying xdata(lags+1:2*lags,:) into 
%      xdata(1:lags,:) are reasonable subsititutes.  These values are used in forming the
%      persistence priors.
[T,nvar]=size(ydata);
nox=isempty(xdata);
if ~nox
   [T2,nx]=size(xdata);
else
   T2=T;nx=0;xdata=zeros(T2,0);
end
% note that x must be same length as y, even though first part of x will not be used.
% This is so that the lags parameter can be changed without reshaping the xdata matrix.
if T2 ~= T, disp('Mismatch of x and y data lengths'),end
X=zeros(T,nvar,lags);
if nargin<4
   nbreaks=0;breaks=[];
else
   nbreaks=length(breaks);
end
breaks=[0;breaks;T];
for nb=1:nbreaks+1
   for i=1:lags
      X(lags+breaks(nb)+1:breaks(nb+1),:,i)=ydata(lags+breaks(nb)+1-i:breaks(nb+1)-i,:);
   end
end
smpl=[];
for nb=1:nbreaks+1
   smpl=[smpl;[breaks(nb)+lags+1:breaks(nb+1)]'];
end
X=[X(:,:) xdata(:,:)];
X=X(smpl,:);
y=ydata(smpl,:);
% Everything now set up with input data for y=Xb+e 
% ------------------Form persistence dummies-------------------
if lambda>0 | mu>0
   ybar=mean(ydata(1:lags,:));
   if ~nox 
      xbar=mean(xdata(1:lags,:));
   else
      xbar=[];
   end
   if lambda>0
      xdum=lambda*[repmat(ybar,1,lags) xbar];
      ydum=zeros(1,nvar);
      ydum(1,:)=lambda*ybar;
      y=[y;ydum];
      X=[X(:,:);xdum];
   end
   if mu>0
      xdum=[repmat(diag(ybar),1,lags) zeros(nvar,nx)]*mu;
      ydum=mu*diag(ybar);
      X=[X;xdum];
      y=[y;ydum];
   end
end
[vl,d,vr]=svd(X(:,:),0);
di=1../diag(d);
B=vl'*y;
B=(vr.*repmat(di',nvar*lags+nx,1))*B;
u=y-X(:,:)*B;
xxi=vr.*repmat(di',nvar*lags+nx,1);
xxi=xxi*xxi';
B=reshape(B,[nvar*lags+nx,nvar]); % rhs variables, equations
By=B(1:nvar*lags,:);
By=reshape(By,nvar,lags,nvar);% variables, lags, equations
By=permute(By,[3,1,2]); %equations, variables, lags to match impulsdt.m
if nox
   Bx=[];
else
   Bx=B(nvar*lags+(1:nx),:)';
end