% MATLAB PROGRAM: coninterval.m

% DESCRIPTION:
% Given an estimated VAR stored in the matrix F, a variance-covariance matrix
% for the residuals (SIGMAV), and a set of estimated correlation of the VAR
% forecast errors at various forecast horizons, this program calculates 
% confidence intervals for the estimted correlations using bootstrapping techniques
% described in Hamilton's Time Series Analysis (pgs. 337-388)

% VARIABLES ASSUMED TO BE IN THE WORKSPACE BEFORE EXECUTION OF comovement.m
% (these will be variables associated with the estimated system)
% ntrend -- number of trends included in the 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
% numlags-- the number of lags included in the estimated 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

% 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

% OUTPUT:
% smplcorrelations-- matrix of correlations for each replication from which the confidence
%                    intervals can be calculated.  A 95% confidence interval for the
%                    estimated correlations can be inferrred from the range of 
%                    correlations in smplcorrelations that includes 95% of the 
%                    values in smplcorrelations.

% 		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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

% Define storage variables for the originally estimated correlations and coefficients
corroutpriceorig=corroutprice;
FORIG=F;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create the simulated data and re-estimate the VAR
% - this procedure is repeated numreplic times
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for nrepl=1:numreplic
   % Create random draws for creating the simulated data
   rdraw=ceil(lvdat*(rand(lvdat,1)));
   % Initialize the simulated data matrix
   simdat=zeros(lvdat,numvars);
   
   % Create time trend variables based on the originally estimated VAR specification
   if ntrend == 1
      trend=(1:1:seriesedate)';
   elseif ntrend == 2
      trend=(1:1:seriesedate)';
      trend2=(trend.*trend);
   else
   end 
   
   % Define the number of right-hand side variables
   numrhs=numlags*numvars+1+ntrend;
   
   % When creating the simulated data, begin by using the first observation
   % and then loop over the observation filling the data matrix of simulated data
   % using the estimated VAR system and random draws from the residuals
   rhs=X(1,:);
   for ii=1:lvdat
      fitval= kron(identity,rhs)*theta;
      simval=fitval'+(RESID(rdraw(ii),:));
      simdat(ii,:)=simval;
      if ii<lvdat
         newrhs=zeros(1,numrhs);
         newrhs(1,1)=1;
         if ntrend==1
            newrhs(1,2)=trend(smplsdate+ii);
            newrhs(1,3:3+numvars-1)=simval;
            newrhs(1,3+numvars:numrhs)=rhs(3:numrhs-numvars);
         elseif ntrend==2
            newrhs(1,2)=trend(smplsdate+ii);
            newrhs(1,3)=trend2(smplsdate+ii);
            newrhs(1,4:4+numvars-1)=simval;
            newrhs(1,4+numvars:numrhs)=rhs(4:numrhs-numvars);
         else
            newrhs(1,2:2+numvars-1)=simval;
            newrhs(1,2+numvars:numrhs)=rhs(2:numrhs-numvars);
         end
         rhs=newrhs;
      else
      end
   end
   
   % Pad the simulated data with out of sample historical data and post sample data
   % using the original data set
   newvardat=zeros(seriesedate,numvars);
   newvardat(1:smplsdate-1,:)=vardat(1:smplsdate-1,:);
   newvardat(smplsdate:smpledate,:)=simdat;
   if smpledate ~= seriesedate
      newvardat(smpledate+1:seriesedate,:)=vardat(smpledate+1:seriesedate,:);
   end
   
   % Estimate the VAR using the new simulated data
   [RESIDN ,thetaN , XN, F, aicN, SIGMAV]=estvar(newvardat,ntrend,numlags);
   % Estimate the correlations of the VAR forecast errors and store the results
   comovement
   smplcorrelations(:,nrepl)=corroutprice;
end