function [A,C,w]=arconf(A,C,w,th)
%ARCONF	Confidence intervals for AR parameter estimates.
%
%	[A,C,w]=ARCONF(A,C,w,th) produces approximate 68.27% 
%	confidence intervals for the estimates of the coefficient 
%	matrix A, the noise covariance matrix C, and the intercept  
%	vector w of an AR(p) model that was fitted using AR. The 
%	input arguments of ARCONF are output of AR.
%
%	The estimated confidence intervals are returned as
%	the imaginary parts of A, C, and w, while the corresponding
%	real parts remain unchanged.
%
%	If no intercept vector w was fitted with AR (i.e., the flag
%	'zero' was an input argument of AR), then [A,C]=ARCONF(A,C,th)
%	produces approximate confidence intervals for the coefficient 
%	matrix A and the noise covariance matrix C only.
%
%	See also: AR

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

m = size(C,1); 				% dimension of state space
p = size(A,2)/m; 			% order of process

gamma = 1; 				% compute 68.27% confidence intervals
sqrt_gamma = sqrt(gamma);

if (nargin == 3)
%  no intercept vector has been fitted
   Aaug = A;
   th = w;
   w = [];
   np = m*p;				% # of parameter vectors of size m
else
   Aaug = [w A];
   np = m*p+1;				% # of parameter vectors of size m
end

dof = th(1,1); 				% # of degrees of freedom in covariance
                                        % estimate

% Get matrix Uinv that appears in the estimates' covariance matrix
Uinv  = th(2:size(th,1),:);
Creal = C;

% Compute confidence intervals for elements of Aaug 
for j=1:m
    for k=1:np
        Aaug(j,k) = Aaug(j,k)+i*sqrt_gamma*sqrt(Uinv(k ,k) * Creal(j,j));
    end
end

if (nargin == 3)
%  No intercept vector has been fitted
   A = Aaug;
else
%  An intercept vector has been fitted => return it and its 
%  confidence intervals separately
   w = Aaug(:,1);
   A = Aaug(:,2:np);
end

% Compute confidence intervals for elements of C
for j=1:m
    for k=1:m
        C(j,k) = Creal(j,k)+i*sqrt_gamma*sqrt(2/dof * Creal(k,k) * Creal(j,j));
    end
end


