function r_k = plot_SACF(et,use_constant_se)
% PLOT_SACF - Plots the sample autocorrelation function with 2*sigma error bars
% 
% Input:
%   et : one-step-ahead error residuals = z_t - \hat{z}_{t-1}(1) 
%   use_constant_se : whether to use a constant standard error for these coefficients O(1/sqrt(n)) 
% 
% Output: 
%   r_k : the sample autocorrelation functions
%   
% Written by:
% -- 
% John L. Weatherwax                2009-04-21
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

if( nargin<2 ) use_constant_se = 0; end

if( exist('/home/weatherw/') ) addpath('/home/weatherw'); end

n     = length(et);
n_r_k = round(0.2*length(et));

r_k    = aasamplebiasedautoc(et,n_r_k); 
r_k(1) = 1.0; % <- override the sample based value with the exact value.

sf=figure; sh_ind=stem( 1:(n_r_k-1), r_k(2:end), 'bx', 'linewidth', 2 ); hold on; % <- we don't explicity plot the (0,r_0=1) pair 
xlabel( 'lag (k)' ); ylabel( 'sample based autocorrelation r_k' ); axis tight; 

% plot the value of the standard errors of the sample autocorrelations: 
if( use_constant_se ) 
  v_r_k  = (1/n) * ones(1,length(r_k)); 
else
  v_r_k  = cumsum( [ 1., 2*r_k(2:end).^2 ] )/n; 
end
se_r_k = sqrt(v_r_k);

figure(sf); hold on; plot( 1:(n_r_k-1), -1.97 * se_r_k(2:end), 'r-' ); 
figure(sf); hold on; plot( 1:(n_r_k-1), +1.97 * se_r_k(2:end), 'r-' ); %axis tight; 
axis( [ 0, n_r_k-0.5, -Inf, +Inf ] ); 

return;