% PROGRAM: calccorr.m

% DESCRIPTION:
% MATLAB program that calculates a concise set of descriptive results about
% the comovement of two variables.  The statistics can be calculated in one of 
% two ways.  The first method calculates the correlations of VAR forecast errors
% for two times-series at various forecast horizons.  The second method calculates
% correlations of filtered time-series at different frequencies using high-pass 
% filters from the frequency domain (proposed by Baxter and King (1994)).
%
% The program requires the following functions and script files:
% 
%		1. pickmodel.m 	(function)
%		2. estvar.m 		(function)
%		3. lags.m 			(function)
%		4. comovement.m 	(script file)
%		5. coninterval.m	(script file)
%		6. highpass.m		(function)
%
% 		These programs and documentation are available at:
%		ftp://weber.ucsd.edu/pub/wdenhaan/comov/
%

% INPUTS:
% See Step 0: USER INPUT For required inputs that must be set prior to execution

% OUTPUTS:
% Method 1: High-Pass Filter
%   corr- correlation of filter series at various frequencies
% Method 2: VAR Forecast Errors
%   corroutprice- correlation of VAR forecast errors at various forecast horizons
%   if confidence intervals are calculated:
%     corroutpriceorig- correlation of VAR forecast errors at various forecast horizons
%								from original data
%		meancorr- average correlation across simulated economies at various forecast horizons
%  	lowerband- lower band of the confidence interval for the specified percentage
%					  confidence interval
%     upperband- upper band of the confidence interval for the specified percentage 
%					  confidence interval
%		sortcorr- correlations for all simulated economies sorted by value (can be used to
%               obtain other percentage confidence intervals
	
% 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

%============================================================================
% STEP 0: USER INPUT
%		-VARIABLES MUST BE SET PRIOR TO EXECUTION
%============================================================================
% Filename of data (this file should contain two columns of data)
infile='c:\matlab\econ220f\internet\testdat.txt';
% What methodology do you wish to use to calculate correlations
% High-pass filter is method 1; VAR forecast errors is method 2
methodology=2;
% Range of forecast horizons for which you wish to view results for VAR method
%  or the periodicity/frequency of interest for High-Pass method 
%  (periodicity should be at least 3 for high-pass filter)

forecastrange=[1 2 3 4 5 6 24 48 84];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% High-Pass method specific user-defined variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Truncation Parameter for the filter
truncparam=40;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% VAR method specific user-defined variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Maximum number of lags to use in selecting best VAR
nlags=12;
% Criterion to use in choosing the best VAR (aic or bic)
infocrit='aic';
% Should a unitroot be imposed in the estimation of the VAR
unitroot=1;
% Should confidence intervals be calculated for the estimated correlations
ciyes=0;
% Confidence interval percentage that should be calculated (a number between 0 and 100)
cipercent=90;
% This number should be divisible by (100-cipercent)/2 
numreplic=2500;

% END STEP 0: USER INPUT
%============================================================================

%============================================================================
% STEP 1: LOAD DATA, CHECK METHODOLOGY, SIZE OF DATA FILE, CONFIDENCE 
%          INTEGER PERCENTAGE
%============================================================================
userdat=load(infile);
if ~(methodology==1 | methodology==2)
   fprintf('%s\n','You have not selected an appropriate methodology, program aborting')     
   break;
end
if size(userdat,2)~=2
   fprintf('%s\n','You do not have exactly two columns in your datafile, program aborting')     
   break;
end
if ciyes==1
   if ~(cipercent>0 & cipercent<100)
      fprintf('%s','Your percentage confidence interval should be between 0 and 100, ')
      fprintf('%s\n','program aborting')     
      break;
   end
   if round(((100-cipercent)/2)*numreplic/100)==0
      fprintf('%s','Your percentage confidence interval and numreplic implies that ')
      fprintf('%s\n','the index for outputting a vector will be zero, program aborting')     
      break;
   end
end
% END STEP 1: CHECK METHODOLOGY
%============================================================================

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% START EXECUTION OF THE PROGRAMS DEPENDING ON METHODOLOGY
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% USE METHOD 1: HIGH-PASS METHOD
if methodology==1
   numvars=2;
   vardat=userdat;
   % Run highpass program
   highpass
   % Print Results
   fprintf('%s\n\n','The correlations you requested are');
   fprintf('%s\n','    Periodicity     Correlation');
   fprintf('%15.4f %15.4f\n',[forecastrange' corr ]');
   
else
   % USE METHOD 2: VAR FORECAST ERRORS METHOD 
   
   %============================================================================
   % STEP 2: CHECK USER INPUTS
   %============================================================================
   if ~(infocrit(1:3)=='aic' | infocrit(1:3)=='bic')
      fprintf('%s','You have not chosen aic or bic for the information criterion, ')
      fprintf('%s\n\n','so the program is aborting')
      break;
   end
   if ~(unitroot==1 | unitroot==0)
      fprintf('%s','You have not chosen a 0 or 1 for whether to impose a unit root, ')
      fprintf('%s\n\n','so the program is aborting')
      break;
   end
   
   % END STEP 2: CHECK USER INPUTS
   %============================================================================
   
   %============================================================================
   % STEP 3: DIFFERENCE DATA IF NECESSARY
   % 		- difference data if unit root is being imposed during estimation
   %		- store data in variable 'vardat'
   %============================================================================
   numvars=2;
   identity=eye(numvars);
   if unitroot==1
      vardat=diff(userdat);
   else
      vardat=userdat;
   end      
   
   %END STEP 3: DIFFERENCE DATA
   %============================================================================
   
   %============================================================================
   % STEP 4: SET SAMPLE RANGE FOR VAR ESTIMATION
   % - define variables seriessdate, seriesedate, smplsdate, smpledate
   % - define the maximum forecast horizon based on user input for which
   %   a correlation will be caluclated
   %============================================================================
   if unitroot==1
      smplsdate=nlags+1;
   else
      smplsdate=nlags+2;
   end
   smpledate=size(vardat,1);
   seriessdate=1;
   seriesedate=size(vardat,1);
   lvdat=smpledate-smplsdate+1;
   datevec=[seriessdate seriesedate smplsdate smpledate];
   maxforchor=max(forecastrange);
   
   % END STEP 4: SET SAMPLE RANGE FOR VAR ESTIMATION
   %============================================================================
   
   %============================================================================
   % STEP 5: PICK THE BEST FITTING MODEL
   % 		- use the procedure 'pickmodel' to select the best model using
   %   	the user selected information criterion
   %============================================================================
   [ntrend, numlags, X, F, RESID, SIGMAV, theta]=pickmodel(nlags,vardat);
   
   % END STEP 5: PICK THE BEST FITTING MODEL
   %=============================================================================
   
   %=============================================================================
   % STEP 6: CALCULATE THE CORRELATIONS 
   % 		- use the script file 'comovement' to calculate the correlations of the
   %   	forecast errors of the two series
   %		- display the results for the user defined forecast horizons
   %=============================================================================
   comovement
   fprintf('%s\n\n','The correlations you requested are');
   fprintf('%s\n','    Forecast Horizon   Correlation');
   fprintf('%15.4f %15.4f\n',[forecastrange' corroutprice(forecastrange) ]');
   
   % END STEP 6: CALCULATE THE CORRELATIONS
   %==============================================================================
   
   %==============================================================================
   % STEP 7: CONFIDENCE INTERVALS
   % 		- use the script file 'coninterval' to calculate the confidence intervals
   %		for the calculated correlation if so requested by user
   %==============================================================================
   if ciyes==1
      smplcorrelations=zeros(maxforchor,numreplic);
      coninterval
      meancorr=mean(smplcorrelations');
      sortcorr=sort(smplcorrelations');
      lowerband=sortcorr(round(((100-cipercent)/2)*numreplic/100),:);
      upperband=sortcorr(round((100-((100-cipercent)/2))*numreplic/100),:);
      fprintf('%s\n\n','The correlations and confidence intervals you requested are');
      fprintf('%s','    Forecast Horizon  Correlation Mean Across Samples ');
      fprintf('%s\n',' Lower Band      Upper Band ');
      fprintf('%15.4f %15.4f %15.4f %15.4f %15.4f\n',[forecastrange'...
            corroutpriceorig(forecastrange) meancorr(forecastrange)'...
            lowerband(forecastrange)' upperband(forecastrange)' ]');
      
   else
   end
   
   % END STEP 7: CONFIDENCE INTERVALS
   %==============================================================================
   
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%==============================================================================
% DOCUMENTATION:
%
% STEPS OF THE ALGORITHM:
% 
% 0. The user defines the following input variables: infile, methodology, 
%		forecastrange, truncparam, nlags, infocrit, unitroot, ciyes, cipercent,
%		and numreplic.  These variables use depends on the methodology selected
%     and may not be used (see STEP 0 for details).
%
% 1. The data specified by infile is loaded and the following user inputs are checked
%    	for useability:
%			- The value of methodology is verified to be either 1==>'High-Pass Method'
%					or 2==>'VAR Forecast Errors Method'
%			- The data file is verified to have 2 columns of data
%		if confidence intervals are to be calculated:
%			- The confidence interval percentage is verified to be between 0 and 100
%
% PROGRAM METHODOLOGY IS SELECTED
% A. High-Pass Method:
%		Data is filtered and correlation of filtered series is calculated for various
%     	frequencies as defined by forecastrange (see program highpass.m for detailed
%			algorithm).  Results are displayed to screen and stored in variable corr.
%
% B. VAR Forecast Method:
%
%	2. User inputs associated with this method are checked:
%			- The value of infocrit is verified to be either 'aic' or 'bic'
%			- The value of unitroot is verified to be either 0==>no or 1==>yes
%
%	3. The data is differenced if unitroot==1 is specified and stored in variable vardat.
%
% 	4. The sample range is defined.  This will depend on the size of infile and the 
%  	maximum number of lags to be used in the VAR estimation.  If a unitroot is 
% 		imposed there will be one less observation lost to differencing the data.  
%		Furthermore, the data will be assumed to start at the same time period for
% 		both the unitroot==0 and unitroot==1 cases.  Therefore the sample start 
%     date (smplsdate) will be the nlags+1 observation for the differenced data and
%		the nlags+2 observation for the level data.  The sample end date (smpledate)
% 		is assumed to be the last observation in infile.  The series date range is also
%     defined and is assumed to be the first and last observations in the datafile and
%  	is stored in the variables seriessdate and seriesedate.
%
% 	5. Using the information criterion specified by the user (either aic or bic) the 
%		best fitting model is chosen.  If a unitroot is not imposed (unitroot==0) then
% 		models including a linear trend and a quadratic trend are considered.  If a 
%		unitroot is imposed (unitroot==1) then models with a linear trend are considered.
% 		In all cases, a constant is included in the estimation.  Also in the estimation
%   	it is imposed that all equations are estimated with the same number of lags.
% 		FUNCTIONS USED IN THIS STEP: pickmodel, estvar, lags  (see these functions for
%     detailed descriptions).
%
% 	6. The correlation between the forecast errors of the two series is calculated.  
%		Using the estimated VAR system, a forecast series is created for each variable 
%		and for each forecast horizon being considered.  The forecast error for each
%		variable is then calculated as the difference between the true and forecasted
%		values and the correlation of these two forecast error series is calculated and
%		reported for the various forecast horizons.
%		SCRIPT FILE comovement.m is used for this step, refer to it for additional
%		documentation.
%
% 	7. If so requested by the user, confidence intervals for the estimated correlation
% 		are calculated using a bootstrapping technique.  Using random draws from the
%  	the estimated VAR residuals (RESID) and the esimated VAR itself, simulated data 
%     is created.	This simulated data is used to estimate a new VAR and calculate a 
%   	new set of correlations.  This is repeated the user-specified number of 
%		times (numreplic) and confidence intervals can be created based on these 
%		results (see Hamilton, Time Series Analysis, 337-338) for more details on the 
% 		algorithm.  
% 		SCRIPT FILE coninterval.m is used for this step
%=======================================================================================
%
% SAMPLE PROGRAM EXECUTION:
% 
%  Set User inputs:
% infile='c:\matlab\internet\testdat.txt';
% methodology=2;
% forecastrange=[1 2 3 4 5 6 24 48 84];
% truncparam=40;
% nlags=12;
% infocrit='aic';
% unitroot=1;
% ciyes=0;
% cipercent=90;
% numreplic=2500;
%
% OUTPUT:
% The minimum aic value is -21.373 found with 12 lags and 
% includes only a constant (ntrend=0) in the estimation of the VAR.
%
% The correlations you requested are
%
%    Forecast Horizon   Correlation
%         1.0000          0.0575
%         2.0000          0.0747
%         3.0000          0.1020
%         4.0000          0.1301
%         5.0000          0.1484
%         6.0000          0.1628
%        24.0000         -0.0544
%        48.0000         -0.2809
%        84.0000         -0.3915
