function y = johansenlags(y_in, lags)
% Use VAR to select lags for Johansen cointegration test


[T,numseries] = size(y_in);
numlags = (1:lags)';
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 = y_in(idxpre,:); % Presample
Y = y_in(idxest,:);  % Estimation sample
for j = 1:numel(numlags)
    Mdl = varm(numseries,j);
%     Mdl.SeriesNames = varnames;
    EstMdl(j) = estimate(Mdl,Y,'Y0',Y0);
    results = summarize(EstMdl(j));
    aic(j) = results.AIC;
end

[~,bestidx] = min(aic);
p = numlags(bestidx);
% BestMdl = EstMdl(bestidx);
% h = gctest(BestMdl,'Display',true);
% 
% res = infer(BestMdl,Y,'Y0',Y0);
% [Response,Lower,Upper] = irf(BestMdl, 'NumObs',60,'Y0',Y0,...
%     'E',res,'Confidence',0.9);

y = p;
end

