% 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