function output = sample_moments(Y,h);
%------------------------------------------------------------
% PURPOSE: calculates the sample multivariate second moments
%          of the time series described by the Txn matrix Y.
%------------------------------------------------------------
% USAGE: output = sample_moments(Y,h)
% where:  Y     = Txn matrix of observations an n-vector y(t)
%         h     = number of desired second moments
%------------------------------------------------------------
% RETURNS: a structure
% output.meth     = 'population moments'
% output.nvar     = no. of variables in y(t)
% output.nobs     = no. of observations on y(t)
% output.nmon     = no. of desired moments
% output.autocov  = 3-tensor of h autocovariances
% output.autocorr = 3-tensor of h autocorrelations
%------------------------------------------------------
% written by: Gregory Givens
%             26 March 2003
%------------------------------------------------------

[T n] = size(Y);
output.meth = 'sample moments';
output.nvar = n; output.nobs = T; output.nmom = h;

AutoCov = zeros(n,n,T);
AutoCorr = zeros(n,n,T);
Ybar = mean(Y);
Y = Y - repmat(Ybar,T,1);
for i = 1:T;
   AutoCov(:,:,i) = (1/(T-i+1))*Y(i:T,:)'*Y(1:T-i+1,:);
end;
D = diag(diag(AutoCov(:,:,1))); Ds = sqrt(D);
for i = 1:T;
    AutoCorr(:,:,i) = inv(Ds)*AutoCov(:,:,i)*inv(Ds);
end;
output.autocov = AutoCov(:,:,1:h+1);
output.autocorr = AutoCorr(:,:,1:h+1);