Merge pull request #5 from lcbarnett/master

v1.3 changes

README.md

@@ -1,7 +1,7 @@
 This is the MVGC Multivariate Granger Causality MATLAB toolbox, previously
 hosted at http://www.sussex.ac.uk/sackler/mvgc
 
-Current version is mvgc_v1.2, last updated Nov. 2021
+Current version is mvgc_v1.3, last updated March 2022
 
 **NOTE**: this version includes an implementation of the new, more efficient and accurate
 state-space method for GC calculation, as detailed in:

core/tsdata_to_var.m

@@ -91,6 +91,20 @@
 [n,m,N] = size(X);
 assert(p < m,'too many lags');
 
+M = N*(m-p); % effective number of observations
+
+if p == 0 % white noise!
+	A = zeros(n,n,0);
+	if nargout > 1
+		X0 = reshape(demean(X),n,M);
+		SIG = (X0*X0')/(M-1);
+		if nargout > 2
+			E = X; % the residuals are just the time series itself!
+		end
+	end
+	return
+end
+
 p1 = p+1;
 pn = p*n;
 p1n = p1*n;
@@ -103,8 +117,6 @@
 
 if  strcmpi(regmode,'OLS') % OLS (QR decomposition)
 
-    M = N*(m-p);
-
     % stack lags
 
     X0 = reshape(X(:,p1:m,:),n,M); % concatenate trials for unlagged observations
@@ -119,8 +131,6 @@
 
     if nargout > 1
         E   = X0-A*XL;             % residuals
-        SIG = (E*E')/(M-1);        % residuals covariance matrix
-        E   = reshape(E,n,m-p,N);  % put residuals back into per-trial form
     end
 
     A = reshape(A,n,n,p);          % so A(:,:,k) is the k-lag coefficients matrix
@@ -177,17 +187,23 @@
 
 	A0 = AF(:,1:n);
     A = reshape(-A0\AF(:,n+1:p1n),n,n,p);
-    if isbad(A); return; end       % something went badly wrong
+    if isbad(A); return; end % something went badly wrong
 
     if nargout > 1
-		M   = N*(m-p);      % residuals lose p lags
-		E   = A0\EF;        % residuals
-        SIG = (E*E')/(M-1); % residuals covariance matrix (unbiased estimator)
-        if nargout > 2      % align residuals per-trial with data (lose p lags)
-			E = cat(2,nan(n,p,N),reshape(E,n,m-p,N));
-		end
+		E = A0\EF; % residuals
     end
 
 else
     error('bad regression mode ''%s''',regmode);
 end
+
+if nargout > 1
+	SIG = (E*E')/(M-1); % residuals covariance matrix (unbiased estimator)
+end
+
+if nargout > 1
+	SIG = (E*E')/(M-1); % residuals covariance matrix (unbiased estimator)
+	if nargout > 2 % pad first p lags of residuals with zeros to align with X
+		E = reshape(E,n,m-p,N);
+	end
+end

startup.m

@@ -13,13 +13,13 @@
 
 global mvgc_version;
 mvgc_version.major = 1;
-mvgc_version.minor = 2;
+mvgc_version.minor = 3;
 
-fprintf('[mvgc startup] Initialising MVGC toolbox version %d.%d\n', mvgc_version.major, mvgc_version.minor);
+fprintf('[MVGC startup] Initialising MVGC toolbox version %d.%d\n', mvgc_version.major, mvgc_version.minor);
 
 %% Set path
 
-% Add mvgc root directory and appropriate subdirectories to path
+% Add MVGC root directory and appropriate subdirectories to path
 
 global mvgc_root;
 mvgc_root = fileparts(mfilename('fullpath')); % directory containing this file
@@ -46,7 +46,7 @@
 % addpath(fullfile(mvgc_root,'testing'));
 % addpath(fullfile(mvgc_root,'maintainer'));
 
-fprintf('[mvgc startup] Added MVGC root directory %s and subdirectories to path\n',mvgc_root);
+fprintf('[MVGC startup] Added MVGC root directory %s and subdirectories to path\n',mvgc_root);
 
 %% Check |mex| files
 
@@ -55,41 +55,41 @@
 global have_genvar_mex;
 have_genvar_mex = exist('genvar_mex','file') == 3;
 if ~have_genvar_mex
-    fprintf(2,'[mvgc startup] WARNING: no ''genvar'' mex file found; please run the ''mvgc_makemex'' script.\n');
-    fprintf(2,'[mvgc startup]          Meanwhile, a slower scripted VAR simulation routine will be used.\n');
+    fprintf(2,'[MVGC startup] WARNING: no ''genvar'' mex file found; please run the ''mvgc_makemex'' script.\n');
+    fprintf(2,'[MVGC startup]          Meanwhile, a slower scripted VAR simulation routine will be used.\n');
 else
-    fprintf('[mvgc startup] All MVGC ''mex'' files for your platform exist\n');
+    fprintf('[MVGC startup] All MVGC ''mex'' files for your platform exist\n');
 end
 
 %% Check for dependencies on other Matlab(R) toolboxes
 
 % Check if we have Statistics toolbox - see if ch2cdf is present
 
 if fexists(@chi2cdf)
-    fprintf('[mvgc startup] Statistics Toolbox(TM) seems to be present.\n');
+    fprintf('[MVGC startup] Statistics Toolbox(TM) seems to be present.\n');
 else
 	addpath(fullfile(mvgc_root,'utils','stats'));
-    fprintf(2,'[mvgc startup] WARNING: Statistics Toolbox(TM) does not seem to be present.\n');
-    fprintf(2,'[mvgc startup]          Will use slower scripted routines (see utils/stats directory).\n');
+    fprintf(2,'[MVGC startup] WARNING: Statistics Toolbox(TM) does not seem to be present.\n');
+    fprintf(2,'[MVGC startup]          Will use slower scripted routines (see utils/stats directory).\n');
 end
 
 % Check if we have Signal Processing toolbox - see if pwelch is present
 
 if fexists(@pwelch)
-    fprintf('[mvgc startup] Signal Processing Toolbox(TM) seems to be present.\n');
+    fprintf('[MVGC startup] Signal Processing Toolbox(TM) seems to be present.\n');
 else
-    fprintf(2,'[mvgc startup] WARNING: Signal Processing Toolbox(TM) does not seem to be present.\n');
-    fprintf(2,'[mvgc startup]          Some spectral estimation routines may not work.\n');
+    fprintf(2,'[MVGC startup] WARNING: Signal Processing Toolbox(TM) does not seem to be present.\n');
+    fprintf(2,'[MVGC startup]          Some spectral estimation routines may not work.\n');
 end
 
 % Check if we have 'dlyap' from the Control System toolbox
 
 if fexists(@dlyap)
-    fprintf('[mvgc startup] Control System Toolbox(TM) seems to be present.\n');
+    fprintf('[MVGC startup] Control System Toolbox(TM) seems to be present.\n');
 else
 	addpath(fullfile(mvgc_root,'utils','control'));
-    fprintf(2,'[mvgc startup] WARNING: Control System Toolbox(TM) does not seem to be present.\n');
-    fprintf(2,'[mvgc startup]          Will use slower scripted routines (see utils/control directory).\n');
+    fprintf(2,'[MVGC startup] WARNING: Control System Toolbox(TM) does not seem to be present.\n');
+    fprintf(2,'[MVGC startup]          Will use slower scripted routines (see utils/control directory).\n');
 end
 
 %% Initialise default random number stream
@@ -100,18 +100,30 @@
 
 rng_seed(-1); % seed from /dev/urandom (Unix/Mac) else from clock
 
-fprintf('[mvgc startup] Random number generator initialised\n');
+fprintf('[MVGC startup] Random number generator initialised\n');
 
 %% Enable all warnings
 
 warning on all
-fprintf('[mvgc startup] All warnings enabled\n');
+fprintf('[MVGC startup] All warnings enabled\n');
+
+% Important notes to users
+
+fprintf('[MVGC startup]\n');
+fprintf('[MVGC startup] NOTE 1: PLEASE DO NOT ADD THE FULL MVGC HIERARCHY TO YOUR MATLAB SEARCH PATH!\n');
+fprintf('[MVGC startup]         Doing so is likely to cause problems. This script has already set up\n');
+fprintf('[MVGC startup]         MVGC paths correctly for your Matlab environment.\n');
+fprintf('[MVGC startup]\n');
+fprintf('[MVGC startup] NOTE 2: It is highly recommended that any single-precision floating-point data\n');
+fprintf('[MVGC startup]         be converted to double precision; some routines may be inaccurate or\n');
+fprintf('[MVGC startup]         numerically unstable for single-precision input.\n');
+fprintf('[MVGC startup]\n');
 
 % Done
 
-fprintf('[mvgc startup] Initialisation complete (you may re-run ''startup'' at any time)\n');
+fprintf('[MVGC startup] Initialisation complete (you may re-run ''startup'' at any time)\n');
 
-fprintf('[mvgc startup] Type ''helpon'' to get started\n');
+fprintf('[MVGC startup] Type ''helpon'' to get started\n');
 
 %%
 % <startup.html back to top>

stats/consistency.m

@@ -58,19 +58,13 @@
 
 X = demean(X);
 
-if N > 1 % multi-trial
-    X = X(:,p+1:m,:);
-    X = X(:,:);              % stack data
-    E = E(:,:);              % stack residuals
-    s = N*(m-p);             % sample size (number of observations)
-else
-    X = X(:,p+1:m);
-    s = m-p;                 % sample size (number of observations)
-end
+M = N*(m-p);                     % effective number of observations
+X = reshape(X(:,p+1:m,:),n,M);   % concatenate trials for data
+E = E(:,:);                      % concatenate trials for residuals
 
-Y = X - E;                   % prediction
+Y = X-E;                         % prediction
 
-Rr = (X*X')/(s-1);           % covariance estimate
-Rs = (Y*Y')/(s-1);           % covariance estimate
+Rr = (X*X')/(M-1);               % covariance estimate
+Rs = (Y*Y')/(M-1);               % covariance estimate
 
 cons = 1 - norm(Rs-Rr)/norm(Rr); % compare matrix norms

stats/rsquared.m

@@ -42,27 +42,22 @@
 function [RSQ,RSQADJ] = rsquared(X,E)
 
 [n,m,N] = size(X);
-p = m-size(E,2);     % number of lags in model
+p = m-size(E,2); % number of lags in model
 assert(m >= n,'too few observations');
 assert(size(E,1) == n && size(E,3) == N,'residuals don''t match data');
 assert(p > 0,'bad number of lags');
 
 X = demean(X);
 
-if N > 1 % multi-trial
-    X = X(:,p+1:m,:);
-    X = X(:,:);              % stack data
-    E = E(:,:);              % stack residuals
-    s = N*(m-p);             % sample size (number of observations)
-else
-    X = X(:,p+1:m);
-    s = m-p;                 % sample size (number of observations)
-end
-r = n*n*p;                   % number of regressors in model
+M = N*(m-p);                   % effective number of observations
+X = reshape(X(:,p+1:m,:),n,M); % concatenate trials for data
+E = E(:,:);                    % concatenate trials for residuals
 
-SSR = sum(E.^2,2)';          % residuals sum of squares
-SST = sum(X.^2,2)';          % total sum of squares
+r = n*n*p;                     % number of regressors in model
+
+SSR = sum(E.^2,2)';            % residuals sum of squares
+SST = sum(X.^2,2)';            % total sum of squares
 
 RSQ = 1 - SSR./SST;
 
-RSQADJ = 1 - ((s-1)/(s-r-1))*(1-RSQ);
+RSQADJ = 1 - ((M-1)/(M-r-1))*(1-RSQ);

stats/whiteness.m

@@ -60,16 +60,21 @@
 function [dw,pval] = whiteness(X,E)
 
 [n,m,N] = size(X);
+p = m-size(E,2);     % number of lags in model
 assert(m >= n,'too few observations');
-assert(size(E,1) == n && size(E,3) == N,'residuals don''t match time series data');
+assert(size(E,1) == n && size(E,3) == N,'residuals don''t match data');
+assert(p > 0,'bad number of lags');
 
 X = demean(X);
 
+M = N*(m-p);                   % effective number of observations
+X = reshape(X(:,p+1:m,:),n,M); % concatenate trials for data
+E = E(:,:);                    % concatenate trials for residuals
+
 dw = zeros(1,n);
 pval = zeros(1,n);
 for i = 1:n
-    Ei = squeeze(E(i,:,:));
-    [dw(i),pval(i)] = durbinwatson(X(:,:),Ei(:)); % concatenate trials (ref. [2])
+    [dw(i),pval(i)] = durbinwatson(X,E(i,:));
 end
 
 function [dw,pval] = durbinwatson(X,E)
@@ -83,7 +88,7 @@
 
 % calculate critical values for the DW statistic using approx method (ref. [1])
 
-A = X*X'; 
+A = X*X';
 B = filter([-1,2,-1],1,X');
 B([1,m],:) = (X(:,[1,m])-X(:,[2,m-1]))';
 D = B/A;