% MATLAB FUNCTION: estvar.m

% USAGE:
% [RESID ,theta , X, F, aic, SIGMAV] = estvar(vardat,numtrend,numlags)

% DESCRIPTION:
% This function uses the following functions:
% lags.m
%
% This procedure is used to estimate a VAR with a specification given by
% the inputs (vardat, numtrend, numlags) and creates the matrix, F, of estimated
% coefficients (VAR(1) structure)

% INPUTS:
% numlags-- the number of lags to include in the VAR estimation
% vardat-- matrix of data (should contain two columns of data) to use in estimation
% numtrend-- a number (0, 1, or 2) indicating whether to include a constant, linear
%            trend or quadratic trend
% numtrend==0 --> VAR has only a constant
% numtrend==1 --> VAR has a constant and linear trend
% numtrend==2 --> VAR has a constant, linear and quadratic trend

% 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:
% RESID --> residuals from VAR estimation
% theta --> the estimated coefficients in a 1-column vector format
% X --> regressor matrix
% F --> estimated coefficient matrix (does not include the contstant
%       or trend coefficients
% aic-- the value of the aic or bic model selection criterion depending on which
%       is being used (used to select best model in the function pickmodel)
% SIGMAV --> variance-covariance matrix of the residuals
%
% 		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 [RESID ,theta , X, F, aic, SIGMAV] = estvar(vardat,numtrend,numlags)

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create time trend variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if numtrend==1
   trend=(1:1:seriesedate)';
elseif numtrend==2
   trend=(1:1:seriesedate)';
   trend2=(trend.*trend);
else
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create matrix of regressors
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X=ones(lvdat,1);
% Fill matrix of regressors
if numtrend==1 
   X=[X trend(smplsdate:smpledate)];
elseif numtrend==2   
   X=[X trend(smplsdate:smpledate) trend2(smplsdate:smpledate)];
end
for i=1:numlags
   X=[X lags(vardat,i,seriessdate, seriesedate,...
         smplsdate, smpledate)];
end

%%%%%%%%%%%%%%%%%%%%%%
% Create regressand
%%%%%%%%%%%%%%%%%%%%%%
Y=vardat(smplsdate:smpledate,:);
Y=Y(:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Estimate Parameters (theta)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
XTEMP=inv(X'*X)*X';
theta=kron(identity,XTEMP)*Y;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create matrices of parameter estimates for each lagged independent variable
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
thetanew=reshape(theta,numvars*numlags+1+numtrend,numvars);
for ii=1:numlags
   eval(['theta' num2str(ii) '=transpose(thetanew(((' num2str(ii)...
         '-1)*numvars+1+numtrend+1:' num2str(ii) '*numvars+1+numtrend),:));']);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create the vector of fitted values using parameter estimates
% Calculate the residuals and then calculate the covariance matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
YFITTED=kron(identity,X)*theta;
RESID=Y-YFITTED;
RESID=reshape(RESID,lvdat,numvars);
SIGMAV=(1/(lvdat))*(RESID'*RESID);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the information criterion
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if infocrit=='bic'
   % BIC criterion should be used
   aic=log(det(SIGMAV))+((log(lvdat)*numvars*(numvars*numlags+numtrend))/lvdat);
else 
   % AIC criterion should be used
   aic=log(det(SIGMAV))+((2*numvars*(numvars*numlags+numtrend))/lvdat);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% If only selecting the best model then stop (onlyaic==1), otherwise continue 
%  to create the state space VAR(1) system of coefficients
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if onlyaic==0
   % Create the State Space VAR(1) notation for the estimated system:
   if unitroot==0
      % Matrix for creating simulated data in coninterval
      F=zeros(numvars*numlags,numvars*numlags);
      for ii=1:numlags
         eval(['F(1:numvars,(' num2str(ii) '-1)*numvars+1:' num2str(ii)...
               '*numvars)=theta' num2str(ii) ';']);
      end
      for iii=1:numlags-1
         F(iii*numvars+1:(iii+1)*numvars,(iii-1)*numvars+1:iii*numvars)=identity;
      end
      
   elseif unitroot==1
      % Reorganize the Estimated VAR system to be in level instead of differences
      % We estimated: Y[t]-Y[t-1]=A1*(Y[t-1]-Y[t-2])+A2*(Y[t-2]-Y[t-3])+...+
      %										AP*(Y[t-p]-Y[t-p-1])
      % This is equivalent to the system:
      %					 Y[t]=(I+A1)*Y[t-1]+(A2-A1)*Y[t-2]+...+(AP-A[P-1])*Y[t-p]+
      %										(-1*AP)*Y[t-p-1]
      
      % Create the State Space VAR(1) notation for the estimated system in levels:
      F=zeros(numvars*(numlags+1),numvars*(numlags+1));
      F(1:numvars,1:numvars)=theta1+identity;
      for ii=2:numlags
         eval(['F(1:numvars,(' num2str(ii) '-1)*numvars+1:' num2str(ii)...
               '*numvars)=theta' num2str(ii) '-theta' num2str(ii-1) ';']);
      end
      eval(['F(1:numvars,numlags*numvars+1:(numlags+1)*numvars)=-1*theta'...
            num2str(numlags) ';']);
      for iii=1:numlags
         F(iii*numvars+1:(iii+1)*numvars,(iii-1)*numvars+1:iii*numvars)=identity;
      end
      
   else
      display('unitroot must be set to 0 or 1');
      break
   end
else 
   F=0;
end