% ----------------------------------------------------------------------- % FUNCTION: regress_window.m % PURPOSE: perform multivariate regression % % INPUT: X - nvar (rows) by nobs (cols) observation matrix % NLAGS - number of lags to include in model % % OUTPUT: % ret.beta - coefficients % ret.u - residuals % ret.RSS - sum-square-error of residuals % ret.Z - covariance matrix of residuals % % Some Original Code Adapted from AKS % Written and Modified by REF Dec 11 2007 % Last Updated December 31 2005 % Ref: Seth, A.K. (2005) Network: Comp. Neural. Sys. 16(1):35-55 % ----------------------------------------------------------------------- function [ret] = regress_window_0305(X,nlags,window_length) % figure regression parameters [nvar,nobs] = size(X); nwindows = nobs / window_length; if(nvar>nobs) error('nvar>nobs, check input matrix'); end if(nlags>window_length) error('nlags>window_length, check parameters'); end %if((nvar*nlags)>(nobs/3)) % warning('Number of observation is not enough, Results may be inaccurate.'); %end for i=1:nvar D = squeeze(reshape(X(i,:),window_length,nwindows)); D = detrend(D); D = D - (ones(window_length,1) * mean(D)); D = D ./ (ones(window_length,1) * std(D)); D = D'; D = D - (ones(nwindows,1) * mean(D)); D = D ./ (ones(nwindows,1) * std(D)); X(i,:) = reshape(D',1,window_length*nwindows); end clear D; % calculate the number of data windows concatonated to input data stream % for example, each row of the data matrix should be a series of seperate data windows which are not necessarily contiguous % [1window 2window .... nwindow] % Now we need to calculate the number of subwindows, each of length nlags (model order) that are in the each window % For example, with a window length of 10 and nlags 4, we would have window_length - nlags subwindows for each window % Original Window [ 1 2 3 4 5 6 7 8 9 10] % % Dependent Variable Values Independent % Subwindow 1 [ 1 2 3 4] 5 % Subwindow 2 [2 3 4 5] 6 % ..... ... % Subwindow 6 [6 7 8 9] 10 % nsubwindows = window_length - nlags; % The number of different case of the model we are going to calculate will be equal to the number of subwindow times the number of windows % construct lag matrices lags = -999*ones([nvar (nsubwindows*nwindows) nlags]); nobs = nsubwindows*nwindows ; parfor ii=1:nlags for kk = 1:(nwindows) lags(1:nvar,(1+((kk-1)*nsubwindows)):(kk*nsubwindows),nlags-ii+1) = X(1:nvar,(ii+((kk-1)*window_length)):(ii+((kk-1)*window_length)+nsubwindows-1)); end end % regression (no constant term) regressors = []; c = zeros(nvar+1,1); for ii=1:nvar regressors = [regressors squeeze(lags(ii,:,:))]; end regressors = squeeze(regressors); XX = regressors' * regressors ; fprintf(1,'\n'); IXX = pinv(XX); IXXX = IXX * regressors' ; clear XX; [a b]=size(regressors); for i=1:a L(i,1)=0.0; end for i=1:a for j=1:b L(i,1)=L(i,1) + (regressors(i,j) * IXXX(j,i)); end end Lcutoff = (2*((nvar*nlags)+1)) / (window_length-nlags+1); FL = ((L - (1/nobs))/(nvar*nlags))./((1-L)/(nobs-(nvar*nlags)-1)) ; alpha=0.05; FLcutoff=finv(1-(alpha/(nobs-(nvar*nlags)-1)),(nvar*nlags), nobs-(nvar*nlags)-1); fprintf(1,'First Loop'); parfor ii=1:nvar xvect = []; for kk=1:(nwindows) xvect = [xvect squeeze(X(ii,((kk-1)*window_length)+(nlags+1):((kk-1)*window_length)+(nsubwindows-1)+(nlags+1)))]; end xvec(ii,:)=xvect ; end parfor ii=1:nvar xdep = xvec(ii,:)'; orig(:,ii) =xdep; beta(:,ii) = IXXX * xdep; pred(:,ii) = regressors*beta(:,ii); u(:,ii) = xdep-regressors*beta(:,ii); RSS1(ii) = sum(u(:,ii).^2); aveerr=RSS1(ii)./nobs; svector(ii,:)=aveerr*diag(IXX); clear aveerr end fprintf(1,'\n\n'); fprintf(1,'Second Loop'); parfor jj=1:nvar eq_inx = setdiff(1:nvar,jj); % vars to include in restricted regression (jj on ii) regressors = []; for kk=1:length(eq_inx) regressors = [regressors squeeze(lags(eq_inx(kk),:,:))]; end regressors = squeeze(regressors); XX = regressors' * regressors ; IXX = pinv(XX); IXXX = IXX * regressors' ; clear IXX; clear XX; for ii=1:nvar if(ii~=jj) xdep = xvec(ii,:)'; beta_r = IXXX * xdep; u_r = xdep-regressors*beta_r; RSS0(ii,jj) = sum(u_r.^2); end end end fprintf(1,'\n\n'); % do Granger f-tests gc = ones(nvar).*NaN; ftest = zeros(nvar); for ii=1:nvar for jj=1:nvar ftest(ii,jj) = ((RSS0(ii,jj)-RSS1(ii))/nlags)/(RSS1(ii)/(nobs-(nlags*nvar)-1)); % causality jj->ii fprob(ii,jj) = 1 - cdff(ftest(ii,jj),nlags,nobs-(nlags*nvar)-1); gc(ii,jj) = log(RSS0(ii,jj)/RSS1(ii)); end end for ii=1:nvar gc(ii,ii) = 0; end for ii=1:nvar-1 for jj=ii+1:nvar UT(:,1) = u(:,ii); UT(:,2) = u(:,jj); M = cov(UT); Rd = det(M); GC_Total(ii,jj) = log((RSS0(ii,jj)*RSS0(jj,ii))/Rd); % causality jj->ii GC_Inter(ii,jj) = log((RSS1(ii)*RSS1(jj))/Rd); % causality jj->ii end end fprintf(1,'finished \n'); % do r-squared df_error = (nobs-1); df_total = (nobs-(nlags*nvar)-1); for ii = 1:nvar xvec = X(ii,nlags+1:end); rss2 = xvec*xvec'; rss(ii) = 1 - (RSS1(ii) ./ rss2); rss_adj(ii) = 1 - ((RSS1(ii)/df_error) / (rss2/df_total) ); end % organize output structure ret.beta = beta; ret.pred = pred; ret.orig = orig; ret.u = u; ret.rss = rss; %ret.rss_adj = rss_adj; ret.Z = cov(u); ret.gc = gc; ret.fs = ftest; ret.fprob = fprob; %ret.GC_Total = GC_Total; %ret.GC_Inter = GC_Inter; %ret.lever = L; %ret.levercutoff = Lcutoff; %ret.fl = FL; %ret.flcutoff = FLcutoff; %ret.svector=svector; %ret.df1=nlags; %ret.df2=nobs-(nlags*nvar)-1; N = nsubwindows*nwindows; nest = nvar*nlags; for ii=1:nvar bic(ii) = -2.00 * log(det(cov(u(ii,:)))) + (log(N) * nest) / N ; aic(ii) = -2.00 * log(det(cov(u(ii,:)))) + (2*nest) / N ; end ret.bic = bic; ret.aic = aic;