function [R, scale]=arqr(v,p,mcor)
%ARQR   QR factorization for AR model.
%
%	[R, SCALE]=ARQR(v,p,mcor) computes the QR factorization needed
%	in the least squares estimation of parameters for an AR(p) 
%	model. If the input flag mcor equals one, a vector of intercept 
%	terms is being fitted. If mcor equals zero, the process v is 
%	assumed to have zero mean. The output argument R is the upper 
%	triangular matrix appearing in the factorization of the AR model, 
%	and SCALE is a vector of scaling factors used to regularize the 
%	factorization. 
%
%	As a lower-level routine, ARQR is called by AR. 

%	Modified 4-Sep-97
%	Author: Tapio Schneider
%		AOS Program, Princeton University 
%		Princeton, NJ 08544-0710
%		tapio@splash.Princeton.EDU

% n: number of time steps; m: dimension of state vectors
[n,m] = size(v);     

ne  = n-p;		% number of block equations of size m
np  = m*p+mcor;		% number of parameter vectors of size m

% Least squares estimation proceeds by solving the LS equations for the 
% linear regression model
%
%                  v(k,:)' = Aaug*u(k,:)' + noise(C)
%
% with Aaug=[w A] and `predictors' 
%
%              u(k,:) = [1 v(k-1,:) ...  v(k-p,:)] 
%
% if the intercept vector w is being fitted. If the process mean is taken 
% to equal zero, then Aaug=A and 
%
%                u(k,:) = [v(k-1,:) ...  v(k-p,:)]. 
%
% The number np therefore is just the dimension of the `predictors' u(k).

% Assemble K matrix (of which QR factorization will be computed)
K = zeros(ne,np+m);	% initialize K
if (mcor == 1)
   % first column of K consists of ones for estimation of intercept vector w
   K(:,1) = ones(ne,1);
end

% Assemble `predictors' u in K 
for j=1:p
    K(:, mcor+m*(j-1)+1:mcor+m*j) = [v(p-j+1:n-j, :)];
end
% Add `observations' v (left hand side of regression model) to K matrix
K(:,np+1:np+m) = [v(p+1:n,:)];

% Compute slightly regularized QR factorization of K
scale = sqrt(eps)*sum(abs(K));
R = triu(qr([K;diag(scale)]));








