%% Adjust Number of Lags for Underlying VAR Model
% By default, gctest assumes an underlying VAR(1) model for all specified response variables.  However, a VAR(1) model might be an inappropriate representation of the data. For example,  the model might not capture all the serial correlation present in the variables.
% To specify a more complex underlying VAR model, you can increase the number of lags by specifying the 'NumLags' name-value pair argument of gctest.
% Consider the Granger causality tests conducted in docid:econ_ug#mw_2703a82b-4f23-4d30-a8db-9a70804b12bf. Load the US macroeconomic data set Data_USEconModel.mat. Convert the M1 money supply and CPI to rates. Convert the US GDP to the real GDP rate.

load Data_USEconModel
m1slrate = price2ret(DataTable.M1SL);
inflation = price2ret(DataTable.CPIAUCSL);
rgdprate = price2ret(DataTable.GDP./DataTable.GDPDEF);

% Preprocess the data by removing all missing observations (indicated by NaN).

idx = sum(isnan([m1slrate inflation rgdprate]),2) < 1;
m1slrate = m1slrate(idx);
inflation = inflation(idx);
rgdprate = rgdprate(idx);
T = numel(m1slrate); % Total sample size

%
% Fit VAR models, with lags ranging from 1 to 4, to the real GDP and inflation rate series. Initialize each fit by specifying the first four observations. Store the Akaike information criteria (AIC) of the fits.

numseries = 2;
numlags = (1:4)';
nummdls = numel(numlags);

% Partition time base.
maxp = max(numlags); % Maximum number of required presample responses
idxpre = 1:maxp;
idxest = (maxp + 1):T;

% Preallocation
EstMdl(nummdls) = varm(numseries,0);
aic = zeros(nummdls,1);

% Fit VAR models to data.
Y0 = [rgdprate(idxpre) inflation(idxpre)]; % Presample
Y = [rgdprate(idxest) inflation(idxest)];  % Estimation sample
for j = 1:numel(numlags)
    Mdl = varm(numseries,numlags(j));
    Mdl.SeriesNames = ["rGDP" "Inflation"];
    EstMdl(j) = estimate(Mdl,Y,'Y0',Y0);
    results = summarize(EstMdl(j));
    aic(j) = results.AIC;
end

p = numlags(aic == min(aic))

% A VAR(3) model yields the best fit.
%
% Assess whether the real GDP rate Granger-causes inflation.  gctest removes p observations from the beginning of the input data to initialize the underlying VAR(p) model for estimation. Prepend only the required p = 3 presample observations to the estimation sample.  Specify the concatenated series as input data. Return the p-value of the test.

rgdprate3 = [Y0((end - p + 1):end,1); Y(:,1)];
inflation3 = [Y0((end - p + 1):end,2); Y(:,2)];
[h,pvalue] = gctest(rgdprate3,inflation3,"NumLags",p)

% The p-value is approximately 0.0008, indicating the existence of strong evidence to reject the null hypothesis of noncausality, that is, that the three real GDP rate lags in the inflation rate equation are jointly zero.  Given the VAR(3) model, there is enough evidence to suggest that the real GDP rate Granger-causes at least one future value of the inflation rate.
%
% Alternatively, you can conduct the same test by passing the estimated VAR(3) model (represented by the varm model object in EstMdl(3)), to the object function gctest. Specify a block-wise test and the "cause" and "effect" series names.

h = gctest(EstMdl(3),'Type',"block-wise",...
    'Cause',"rGDP",'Effect',"Inflation")

% Copyright 2018 The MathWorks, Inc.