function [pval,F,sig,pvalun,Fun,sigun,A,morder] = fmvgc(X,alpha,morder,Ae)

%% fmvgc
%
% p-values for sample MVGC based on theoretical asymptotic null distribution
%
% <matlab:open('fmvgc.m') code>
%
%% Syntax
%
%     [pval,F,sig,pvalun,Fun,sigun,A,morder] = fmvgc(X,alpha,morder,Ae)
%
%% Arguments
%
% See also <mvgchelp.html#4 Common variable names and data structures>.
%
% _input_
%
%     X          time series data
%     alpha      significance level for null-hypothesis testing
%     morder     model order to use ('actual', 'AIC', 'BIC' or supplied numerical value)
%     Ae         [optional] externally estimated AR model.
%
% _output_
%
%     pval       matrix of p-values
%     F          Granger causality matrix
%     sig        significance of causality with correction
%     pvalun     matrix of unconditional MVGC p-values
%     Fun        unconditional Granger causality matrix
%     sigun      significance of unconditional causality w correction
%     A          autocorrelation matrix estimate
%     morder     actual model order
% Parameters
%%

ntrials   = 1;      % number of trials
%nobs      = ?;     % number of observations per trial

regmode   = 'LWR';  % VAR model estimation regression mode ('OLS', 'LWR' or empty for default)
icregmode = 'LWR';  % information criteria regression mode ('OLS', 'LWR' or empty for default)
if nargin < 3,
    morder    = 'AIC';  % model order to use ('actual', 'AIC', 'BIC' or supplied numerical value)
end;
momax     = 10;     % maximum model order for model order estimation

%acmaxlags = 1000;  % maximum autocovariance lags 
acmaxlags = [];      % (empty for automatic calculation)

tstat     = [];     % statistical test for MVGC:  'F' for Granger's F-test (default) or 'chi2' for Geweke's chi2 test

%alpha     = 0.01;   % significance level for significance test

mhtc      = 'FDR';  % multiple hypothesis test correction (see routine 'significance')

fs        = 1;    % sample rate (Hz)
fres      = [];     % frequency resolution (empty for automatic calculation)

%% VAR data (A3)

[nvars,nobs]=size(X); % Each variable row-wise oriented

if nvars > nobs,
    X=X';
    [nvars,nobs]=size(X);
end;


% Model order estimation (A2)
% Calculate information criteria up to specified maximum model order.

[AIC,BIC,moAIC,moBIC] = tsdata_to_infocrit(X,momax,icregmode);

fprintf('\nbest model order (AIC) = %d\n',moAIC);
fprintf('best model order (BIC) = %d\n',moBIC);

if nargin > 3,
    amo = size(Ae,3); % actual model order
    fprintf('actual model order     = %d\n',amo);
    
    if  strcmpi(morder,'actual')
        morder = amo;
        fprintf('\nusing actual model order = %d\n',morder);
    end;
else
    if  strcmpi(morder,'actual'),
        morder = 'AIC';
        fprintf('\nNo actual model info, ');
    end;
end;


% Select model order.

if  strcmpi(morder,'actual')
    morder = amo;
    fprintf('\nusing actual model order = %d\n',morder);
elseif strcmpi(morder,'AIC')
    morder = moAIC;
    fprintf('\nusing AIC best model order = %d\n',morder);
elseif strcmpi(morder,'BIC')
    morder = moBIC;
    fprintf('\nusing BIC best model order = %d\n',morder);
else
    fprintf('\nusing specified model order = %d\n',morder);
end

% % VAR model estimation (A2)

% Estimate VAR model of selected order from data.

[A,SIG] = tsdata_to_var(X,morder,regmode);

% Check for failed regression

assert(~isbad(A),'VAR estimation failed');

%If provided, externally estimated A will be used.
if nargin > 3,   A=Ae;  end;

% NOTE: at this point we have a model and are finished with the data! - all
% subsequent calculations work from the estimated VAR parameters A and SIG.

% % Autocovariance calculation (A5)

% The autocovariance sequence drives many Granger causality calculations (see
% next section). Now we calculate the autocovariance sequence G according to the
% VAR model, to as many lags as it takes to decay to below the numerical
% tolerance level, or to acmaxlags lags if specified (i.e. non-empty).

[G,info] = var_to_autocov(A,SIG,acmaxlags);

% The above routine does a LOT of error checking and issues useful diagnostics.
% If there are problems with your data (e.g. non-stationarity, colinearity,
% etc.) there's a good chance it'll show up at this point - and the diagnostics
% may supply useful information as to what went wrong. It is thus essential to
% report and check for errors here.

var_info(info,true); % report results (and bail out on error)

% % Granger causality calculation: time domain (A13)

% Calculate time-domain pairwise-conditional causalities - this just
% requires the autocovariance sequence.
F = autocov_to_pwcgc(G);

% Check for failed GC calculation
assert(~isbad(F,false),'GC calculation failed');

% Significance test using theoretical null distribution, 

pval = mvgc_pval(F,morder,nobs,ntrials,1,1,nvars-2,tstat);% take careful
%                                                         % note of arguments!

% and adjusting for multiple hypotheses.

sig  = significance(pval,alpha,mhtc);

%%
% Unconditional multivariate Granger causality

[G2,info2] = var_to_autocov(A,SIG,acmaxlags); %,acdectol,aitr,maxiters,maxrelerr)

for x = 1:nvars,
    for y = 1:nvars,
        
        if x ~= y
            F2 = autocov_to_mvgc(G2,x,y);
            Fun(x,y) = F2;
        else
            Fun(x,y)=NaN;
        end;
    end;
end;

nx = 1; ny = 1; nz = 0; N = 1;

pvalun = mvgc_pval(Fun,morder,nobs,N,nx,ny,nz); %,tstat)
sigun  = significance(pvalun,alpha,mhtc);

%%
% % From MVGC Toolbox v1.0. by Lionel Barnett and Anil K. Seth, 2012.
% % See file license.txt in installation directory for licensing terms.