Одного кирпичика для грааля не хватает (вопрос к квантам)

Нашел текст теста на собачьем языке R:

Granger causality test

Model 1: fii ~ Lags(fii, 1:5) + Lags(rM, 1:5)
Model 2: fii ~ Lags(fii, 1:5)
Res.Df Df F Pr(>F)
1 629
2 634 5 2.5115 0.02896 *
— Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Granger causality test

Model 1: rM ~ Lags(rM, 1:5) + Lags(fii, 1:5)
Model 2: rM ~ Lags(rM, 1:5)
Res.Df Df F Pr(>F)
1 629
2 634 5 1.1804 0.3172

Надо либо с этим R разбираться, есть так компилятор и ли че… Или на нормальный переписать. C# или C++ или хоть на Бейсик.

нашелся код для матлаба (хоть смысл могу понять):
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: [email protected]
% $Revision: 1.0.0 $ $Date: 09/30/2009 $
% $Revision: 1.0.1 $ $Date: 10/20/2009 $

% References:
% [1] Granger, C.W.J., 1969. «Investigating causal relations by econometric
% models and cross-spectral methods». Econometrica 37 (3), 424?438.

%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

x_lag = find(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

y_lag = find(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)));