% MATLAB FUNCTION: pickmodel.m

% USAGE:
% [ntrend, minlag, X, F, RESID, SIGMAV, theta] = pickmodel(nlags,vardat)

% DESCRIPTION:
% This function uses the following functions:
% estvar.m
%
% This function estimates various VAR models (using the function estvar) and 
% then chooses the best performing model based on either the aic or bic 
% information criterion.
%
% Estimated VARs include:
% 1. only a constant
% 2. a constant and linear trend
% 3. a constant, linear trend, and quadratic trend (unitroot==0 case only)
%
% For each of the possible VAR cases above, a VAR is estimated with 1 lag,
% 2 lags, ... up to nlags lags.  The program imposes that each equation have
% the same number of lags.

% INPUTS:
% nlags-- the maximum number of lags to allow for in VAR estimation
% vardat-- matrix of data (should contain two columns of data)

% GLOBAL VARIABLES ASSUMED TO BE DEFINED IN WORKSPACE:
% datevec- vector that gives series and sample start and end indices
% numvars-- number of variables (should be set=2)
% seriessdate-- index of the series start date (should = 1)
% seriesedate-- index of the series end date (should be the index of the last element
%               in the time-series)
% smplsdate-- index of the first element in the sample
% smpledate-- index of the last element in the sample
% identity-- a numvars by numvars indentity matrix
% lvdat-- number of observations in the sample range (smpledate-smplsdate+1)
% infocrit-- a string value equal to either 'aic' or 'bic'; gives the type of 
%            model selection criterion to use in selecting best VAR model
% unitroot-- value of 0 or 1, indicates whether a unitroot should be imposed in estimation
% maxforchor-- maximum periods ahead to estimate a correlation (maximum of forecastrange)
% numreplic-- number of replicated economies to use in confidence interval calculation
%             (used in coninterval.m program)
% onlyaic-- value of 0 or 1, indicates whether only selecting best model or actually 
%           estimating coefficients (allows program to skip steps if only selecting
%           best model

% OUTPUTS:
% The function returns the following values:
% ntrend --> number of trends included in best fitting VAR.
% ntrend==0 --> VAR has only a constant
% ntrend==1 --> VAR has a constant and linear trend
% ntrend==2 --> VAR has a constant, linear and quadratic trend
% minlag --> the number of lags in the best fitting VAR
% X --> regressor matrix
% F --> estimated coefficient matrix (does not include the contstant
%       or trend coefficients
% RESID --> residuals from VAR estimation
% SIGMAV --> variance-covariance matrix of the residuals
% theta --> the estimated coefficients in a 1-column vector format
%
% 		These programs are available at:
%		ftp://weber.ucsd.edu/pub/wdenhaan/comov/
%
% September 25, 2000	
% Written by Steve Sumner, University of California, San Diego
% Support provided by Wouter den Haan's NSF grant #9708587 is gratefully acknowledged.
%
% Documentation of the program is given at the end of the program
% and on the website http://weber.ucsd.edu/~wdenhaan/soft.html
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ntrend, minlag, X, F, RESID, SIGMAV, theta] = pickmodel(nlags,vardat)

global datevec numvars seriessdate seriesedate smplsdate smpledate identity...
   lvdat infocrit unitroot maxforchor numreplic onlyaic

onlyaic=1; % if equal to 1, then only calculate information criterion

%##########################################################################
% Loop over the maximum possible number of lags to include VAR estimation
%##########################################################################
for nlagloop=1:nlags
   
   numlags=nlagloop;  	% number of lags to include in VAR estimation
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % ESTIMATE VAR WITH NO TRENDS (numtrend==0)
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   numtrend=0;
   [RESID ,theta , X, F, aic, SIGMAV]=estvar(vardat,numtrend,numlags);
   
   % Check to see if this VAR specification has minimized the information criterion
   %  if it does, store the result and keep track of this specification
   if nlagloop==1
      minaic=aic;
      minlag=nlagloop;
      ntrend=0;   
   end
   if aic<minaic
      minaic=aic;
      minlag=nlagloop;
      ntrend=0;   
   end
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % ESTIMATE VAR WITH A LINEAR TREND (numtrend==1)
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   numtrend=1;
   [RESID ,theta , X, F, aic, SIGMAV]=estvar(vardat,numtrend,numlags);
   
   % Check to see if this VAR specification has minimized the information criterion
   %  if it does, store the result and keep track of this specification
   if aic<minaic
      minaic=aic;
      minlag=nlagloop;
      ntrend=1;   
   end
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % ESTIMATE VAR WITH A LINEAR AND QUADRATIC TREND (numtrend==2)
   % **ONLY CONSIDER THIS SPECIFICATION IF A UNITROOT HAS NOT BEEN IMPOSED
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   if unitroot==0
      numtrend=2;
      [RESID ,theta , X, F, aic, SIGMAV]=estvar(vardat,numtrend,numlags);
      
      % Check to see if this VAR specification has minimized the information criterion
      %  if it does, store the result and keep track of this specification
      if aic<minaic
         minaic=aic;
         minlag=nlagloop;
         ntrend=2;
      end
   end
end %nlag loop

% Define text variables for output of the winning specification
if ntrend==0
   ntrendword='only a constant';
elseif ntrend==1
   ntrendword='a linear trend';
else
   ntrendword='a linear and quadratic trend';
end

% Print out results of the model specification search
text1=['The minimum aic value is ' num2str(minaic) ' found with ' num2str(minlag) ' lags and '];
text2=['includes ' ntrendword ' (ntrend=' num2str(ntrend) ') in the estimation of the VAR.'];
fprintf(1,'\n\n%s\n',text1);
fprintf(1,'%s\n\n',text2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ESTIMATE THE 'BEST' MODEL THAT WAS SELECTED ABOVE
% - return the values of RESID,theta,lvdat,X,F,aic to 
%   the calling program frontend.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
onlyaic=0;
[RESID ,theta , X, F, aic, SIGMAV]=estvar(vardat,ntrend,minlag);
