function [fitted, b, ts, r2, rss, BIC, df]=ols_reg_arx(y, x, lag_y, lag_x, c)
%OLSH(X, Y)      	runs a linear regression of the vector y on
%                  	the matrix x and reports the results.
%c=1 includes the constant

%COMPUTE NEW MATRIX OF LAGS Z%
Z_x = [];
for j = 1:lag_x;
    Z_x = [Z_x x(lag_x+1-j:end-j,1) ];
end;

Z_y = [];

for j = 1:lag_y;
    Z_y = [Z_y y(lag_y+1-j:end-j,1) ];
end;

mm=max([lag_x lag_y]);

Z=[Z_y(mm-lag_y+1:end,:) Z_x(mm-lag_x+1:end,:)];
Y=y(mm+1:end,1);


if c == 1
    Z= [ones(size(y(mm+1:end,1),1),1) Z];
else
    Z=Z;
end

%COMPUTE LENGTH OF NEW THE VECTOR AND MATRIX%
[t, k] = size(Z);        % t = # of obs., k = # of regressors
[m, n] = size(Y);
df = t - k;              % degrees of freedom


%POSSIBLE ERROR MESSAGES%
if n > 1
    error('y must be a vector.');
end
if size(y,1) ~= size(x,1)
    error('y and x must have the same number of rows.');
end
if k > size(y(mm+1:end,1),1);
    error('More parameters than observations.');
end

if lag_y==0 & lag_x ==0 & c==0
    error('no regressors.');
end


%--- COMPUTE OLS STATISTICS ---%
b    = inv(Z'*Z)*Z'*Y;                       % coefficients
u    = Y - Z*b;                              % residuals
tss  = (t-1) * std(Y)^2;                     % total sum of squares
rss   = u'*u;                                % residual sum of squares
ess   = tss - rss;                           % explained sum of squares


%--- compute the variance/covariance matrix of b ---%
vc = rss * inv(Z'*Z) / df;

%--- compute various regression statistics ---%
se = sqrt(diag(vc));                       % standard errors of b
ts = b ./ se;                              % t statistics
r2 = ess / tss;                            % R-squared
dw = sum( (u(2:t) - u(1:t-1)).^2 )' / rss; % Durbin-Watson statistic

BIC =log(rss/t)+k*log(t)/t; 
fitted = Z*b;