Список проголосовавших за комментарий "нашелся код для матлаба (хоть смысл могу понять):<br/> function [F,c_v] = granger_cause(x,y,alpha,max_lag)<br/> % [F,c_v] = granger_cause(x,y,alpha,max_lag)<br/> % Granger Causality test<br/> % Does Y Granger Cause X?<br/> %<br/> % User-Specified Inputs:<br/> % x — A column vector of data<br/> % y — A column vector of data<br/> % alpha — the significance level specified by the user<br/> % max_lag — the maximum number of lags to be considered<br/> % User-requested Output:<br/> % F — The value of the F-statistic<br/> % c_v — The critical value from the F-distribution<br/> %<br/> % The lag length selection is chosen using the Bayesian information<br/> % Criterion <br/> % Note that if F > c_v we reject the null hypothesis that y does not<br/> % Granger Cause x<br/> <br/> % Chandler Lutz, UCR 2009<br/> % Questions/Comments: chandler.lutz@email.ucr.edu<br/> % $Revision: 1.0.0 $ $Date: 09/30/2009 $<br/> % $Revision: 1.0.1 $ $Date: 10/20/2009 $<br/> <br/> % References:<br/> % [1] Granger, C.W.J., 1969. «Investigating causal relations by econometric<br/> % models and cross-spectral methods». Econometrica 37 (3), 424?438.<br/> <br/> %Make sure x & y are the same length<br/> if (length(x) ~= length(y))<br/> error('x and y must be the same length');<br/> end<br/> <br/> %Make sure x is a column vector<br/> [a,b] = size(x);<br/> if (b>a)<br/> %x is a row vector — fix this<br/> x = x';<br/> end<br/> <br/> %Make sure y is a column vector<br/> [a,b] = size(y);<br/> if (b>a)<br/> %y is a row vector — fix this<br/> y = y';<br/> end<br/> <br/> %Make sure max_lag is >= 1<br/> if max_lag < 1<br/> error('max_lag must be greater than or equal to one');<br/> end<br/> <br/> %First find the proper model specification using the Bayesian Information<br/> %Criterion for the number of lags of x<br/> <br/> T = length(x);<br/> <br/> BIC = zeros(max_lag,1);<br/> <br/> %Specify a matrix for the restricted RSS<br/> RSS_R = zeros(max_lag,1);<br/> <br/> i = 1;<br/> while i <= max_lag<br/> ystar = x(i+1:T,:);<br/> xstar = [ones(T-i,1) zeros(T-i,i)];<br/> %Populate the xstar matrix with the corresponding vectors of lags<br/> j = 1;<br/> while j <= i<br/> xstar(:,j+1) = x(i+1-j:T-j);<br/> j = j+1;<br/> end<br/> %Apply the regress function. b = betahat, bint corresponds to the 95%<br/> %confidence intervals for the regression coefficients and r = residuals<br/> [b,bint,r] = regress(ystar,xstar);<br/> %Find the bayesian information criterion<br/> BIC(i,:) = T*log(r'*r/T) + (i+1)*log(T);<br/> %Put the restricted residual sum of squares in the RSS_R vector<br/> RSS_R(i,:) = r'*r;<br/> i = i+1;<br/> end<br/> <br/> x_lag = find(min(BIC));<br/> <br/> %First find the proper model specification using the Bayesian Information<br/> %Criterion for the number of lags of y<br/> <br/> BIC = zeros(max_lag,1);<br/> <br/> %Specify a matrix for the unrestricted RSS<br/> RSS_U = zeros(max_lag,1);<br/> <br/> i = 1;<br/> while i <= max_lag<br/> ystar = x(i+x_lag+1:T,:);<br/> xstar = [ones(T-(i+x_lag),1) zeros(T-(i+x_lag),x_lag+i)];<br/> %Populate the xstar matrix with the corresponding vectors of lags of x<br/> j = 1;<br/> while j <= x_lag<br/> xstar(:,j+1) = x(i+x_lag+1-j:T-j,:);<br/> j = j+1;<br/> end<br/> %Populate the xstar matrix with the corresponding vectors of lags of y<br/> j = 1;<br/> while j <= i<br/> xstar(:,x_lag+j+1) = y(i+x_lag+1-j:T-j,:);<br/> j = j+1;<br/> end<br/> %Apply the regress function. b = betahat, bint corresponds to the 95%<br/> %confidence intervals for the regression coefficients and r = residuals<br/> [b,bint,r] = regress(ystar,xstar);<br/> %Find the bayesian information criterion<br/> BIC(i,:) = T*log(r'*r/T) + (i+1)*log(T);<br/> RSS_U(i,:) = r'*r;<br/> i = i+1;<br/> end<br/> <br/> y_lag = find(min(BIC));<br/> <br/> %The numerator of the F-statistic<br/> F_num = ((RSS_R(x_lag,:) — RSS_U(y_lag,:))/y_lag);<br/> <br/> %The denominator of the F-statistic<br/> F_den = RSS_U(y_lag,:)/(T-(x_lag+y_lag+1));<br/> <br/> %The F-Statistic<br/> F = F_num/F_den;<br/> <br/> c_v = finv(1-alpha,y_lag,(T-(x_lag+y