% GRANGER CAUSALITY NONSTANDARD DISTRIBUTION WHEN BOTH PROCESSES ARE 
% DRIFTLESS I(1)
% This code uses 
%--------------------------------------------------------------------------
% Note 1: Requires JPLeSAge toolbox [http://www.spatial-econometrics.com/]

function sim_theorem1
m=10000;
T=10000;
K=10000;
for i=1:m
Sim=(100*i/m)
    F_GRANGER(i)=GC_SD(1,1,T);
    FG_ASINT(i)=FG_TEO(K);
clc;
end
[a,b]=ksdensity(F_GRANGER);
[c,d]=ksdensity(FG_ASINT);
figure(1)
subplot(211)
 plot(b,a),hold on, plot(d,c,'r');
 subplot(212)
 cdfplot(F_GRANGER),hold on, cdfplot(FG_ASINT)
 

function F_GRANGER=GC_SD(sigmax,sigmay,T)

 Ux=normrnd(0,sigmax,T,1);
 Uy=normrnd(0,sigmay,T,1);

 y=cumsum(Uy);
 x=cumsum(Ux);

y_lag=y(1:T-1);
y_con=y(2:T);

x_lag=x(1:T-1);
x_con=x(2:T);

results=ols(y_con,[y_lag,x_lag]);
rcuad_aux= results.resid;
rcuad=rcuad_aux.^2;

results_r=ols(y_con,[y_lag]);
rcuad_raux= results_r.resid;
rcuad_r=rcuad_raux.^2;

RRSS=sum(rcuad_r);
URSS=sum(rcuad);

F_GRANGER= (RRSS-URSS)/(URSS/(T-1));


function [FG_ASINT]=FG_TEO(T)

y=randn(T,1);Ey=cumsum(y);
x=randn(T,1);Ex=cumsum(x);

E2y=sum(Ey.^2)/(T^2);
E2x=sum(Ex.^2)/(T^2);
EUy=sum(Ey.*y)/(T);
Eyx=sum(Ey.*Ex)/(T^2);
ExUy=sum(Ex.*y)/(T);
U2y=sum(y.^2)/(T);
FG_ASINT=(((E2y*ExUy)-(EUy*Eyx))^2)/(E2y*U2y*(E2x*E2y-(Eyx^2)));
