xL(y, x, lags)

function [xx, xy] = xL(y, x, lags)
%xL Generates X'*x and X'*y for ARX(n) regression
%   [xx, xy] = xL(y, x, lags) returns X'*X and X'*y matrices for the
%   regression y_t = \sum_{i=1}^N \sum_{k=1}^n_i \beta_{ik} x_{i, t-k}
%   where:
%       x_i - columns of the input matrix 'x'
%       n_i - ith element of the input vector 'lags'
%       \beta_{ik} - coefficients to estimate
%
%   This is the version specifically optimized for very large regressions,
%   with number of observations >1e5.
%
%example:
%       run regression of y on 200 lags of y and 100 lags of y.^2
%       [xtx, xty] = xL(y(2:end), [y(1:end-1) y(1:end-1).^2], [200 100]);
%       beta = inv(xtx)*xty;
%
%   see also REGRESS, LAGMATRIX (from UCSD GARCH toolbox)
%   For bug reports please contact [e-mail] 

    assert(size(x, 2) == length(lags));
    assert(size(y, 1) == size(x, 1));
    
    N = cols(x);
    xx = zeros(sum(lags));
    xy = zeros(sum(lags), 1);
    
    for i = 1:N
        for j = i:N
            offseti = sum(lags(1:i-1));
            offsetj = sum(lags(1:j-1));
            if (j == i)
                [t1, xy(offseti+1:offseti+lags(i)) ] = auc_xtx_symm(x(:,i), lags(i), y);
                xx(offseti+1:offseti+lags(i), offsetj+1:offsetj+lags(j)) = toeplitz(t1);
            else
                [t1, t2] = auc_xtx(x(:, i), lags(i), x(:, j), lags(j));
                warning off all
                temp = toeplitz(t1, t2);
                warning on all
                xx(offseti+1:offseti+lags(i), offsetj+1:offsetj+lags(j)) = temp;
                xx(offsetj+1:offsetj+lags(j), offseti+1:offseti+lags(i)) = temp';
            end
        end
    end
end