function [F,c_v] = granger_cause(x,y,alpha,max_lag)

% [F,c_v] = granger_cause(x,y,alpha,max_lag)

% Granger Causality test

% Does Y Granger Cause X?

%

% User-Specified Inputs:

%   x -- A column vector of data

%   y -- A column vector of data

%   alpha -- the significance level specified by the user

%   max_lag -- the maximum number of lags to be considered

% User-requested Output:

%   F -- The value of the F-statistic

%   c_v -- The critical value from the F-distribution

%

% The lag length selection is chosen using the Bayesian information

% Criterion 

% Note that if F > c_v we reject the null hypothesis that y does not

% Granger Cause x

% Chandler Lutz, UCR 2009

% Questions/Comments: chandler.lutz@email.ucr.edu

% $Revision: 1.0.0 $  $Date: 09/30/2009 $

% $Revision: 1.0.1 $  $Date: 10/20/2009 $

% $Revision: 1.0.2 $  $Date: 03/18/2009 $

% References:

% [1] Granger, C.W.J., 1969. "Investigating causal relations by econometric

%     models and cross-spectral methods". Econometrica 37 (3), 424–438.

% Acknowledgements:

%   I would like to thank Mads Dyrholm for his helpful comments and

%   suggestions

%Make sure x & y are the same length

if (length(x) ~= length(y))

    error('x and y must be the same length');

end

%Make sure x is a column vector

[a,b] = size(x);

if (b>a)

    %x is a row vector -- fix this

    x = x';

end

%Make sure y is a column vector

[a,b] = size(y);

if (b>a)

    %y is a row vector -- fix this

    y = y';

end

%Make sure max_lag is >= 1

if max_lag < 1

    error('max_lag must be greater than or equal to one');

end

%First find the proper model specification using the Bayesian Information

%Criterion for the number of lags of x

T = length(x);

BIC = zeros(max_lag,1);

%Specify a matrix for the restricted RSS

RSS_R = zeros(max_lag,1);

i = 1;

while i <= max_lag

    ystar = x(i+1:T,:);

    xstar = [ones(T-i,1) zeros(T-i,i)];

    %Populate the xstar matrix with the corresponding vectors of lags

    j = 1;

    while j <= i

        xstar(:,j+1) = x(i+1-j:T-j);

        j = j+1;

    end

    %Apply the regress function. b = betahat, bint corresponds to the 95%

    %confidence intervals for the regression coefficients and r = residuals

    [b,bint,r] = regress(ystar,xstar);

    %Find the bayesian information criterion

    BIC(i,:) = T*log(r'*r/T) + (i+1)*log(T);

    %Put the restricted residual sum of squares in the RSS_R vector

    RSS_R(i,:) = r'*r;

    i = i+1;

end

[dummy,x_lag] = min(BIC);

%First find the proper model specification using the Bayesian Information

%Criterion for the number of lags of y

BIC = zeros(max_lag,1);

%Specify a matrix for the unrestricted RSS

RSS_U = zeros(max_lag,1);

i = 1;

while i <= max_lag

    ystar = x(i+x_lag+1:T,:);

    xstar = [ones(T-(i+x_lag),1) zeros(T-(i+x_lag),x_lag+i)];

    %Populate the xstar matrix with the corresponding vectors of lags of x

    j = 1;

    while j <= x_lag

        xstar(:,j+1) = x(i+x_lag+1-j:T-j,:);

        j = j+1;

    end

    %Populate the xstar matrix with the corresponding vectors of lags of y

    j = 1;

    while j <= i

        xstar(:,x_lag+j+1) = y(i+x_lag+1-j:T-j,:);

        j = j+1;

    end

    %Apply the regress function. b = betahat, bint corresponds to the 95%

    %confidence intervals for the regression coefficients and r = residuals

    [b,bint,r] = regress(ystar,xstar);

    %Find the bayesian information criterion

    BIC(i,:) = T*log(r'*r/T) + (i+1)*log(T);

    RSS_U(i,:) = r'*r;

    i = i+1;

end

[dummy,y_lag] =min(BIC);

%The numerator of the F-statistic

F_num = ((RSS_R(x_lag,:) - RSS_U(y_lag,:))/y_lag);

%The denominator of the F-statistic

F_den = RSS_U(y_lag,:)/(T-(x_lag+y_lag+1));

%The F-Statistic

F = F_num/F_den;

c_v = finv(1-alpha,y_lag,(T-(x_lag+y_lag+1)));