% PROGRAM: highpass.m
%
% MATLAB program that filters two series for user-specified periodicities and
% then calculates the correlation between the two series at each of the 
% periodicities.

% REQUIRED INPUTS (variables that must be defined before execution):
% forecastrange-- periodicities of interest (measured in months, quarters, etc.)
% truncparam-- truncation parameter for filtering (truncated at lead/lag K)
% vardat-- matrix of data (needs to be two columns, one column for each time-series)
% numvars-- number of variables in system (should be set=2)

% September 25, 2000	
% Written by Steve Sumner, University of California, San Diego
% Support provided by Wouter den Haan's NSF grant is gratefully acknowledged.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Initialize the correlation storage vector
corr=zeros(size(forecastrange,2),1);

% Loop over the desired periodicities and calculate the correlations
for nn=1:size(forecastrange,2)
   
   % Initialize the weighting coefficients vectors
   b=zeros(2*truncparam+1,1);
   a=zeros(2*truncparam+1,1);
   
   % Define the frequency based on the periodicities
   w=2*pi/forecastrange(nn);
   
   % Define the weighting matrix
   for i=1:truncparam
      b(truncparam+1-i)=-sin(w*i)/(pi*i);
      b(i+truncparam+1)=-sin(w*-i)/(pi*-i);
   end
   b(truncparam+1)=1-w/pi;
   
   % Adjust the weighting matrix due to truncation
   theta=-1*sum(b)/(2*truncparam+1);
   a=b+theta;
   
   % Pull off the individual unfiltered times-series
   series1=vardat(:,1);
   series2=vardat(:,2);
   
   % Define the length of the original time-series (# of observations)
   lvdat=(size(series1,1));
   
   % Initialize the storage vectors for the filtered series
   filterseries1=zeros(lvdat-2*truncparam,1);
   filterseries2=zeros(lvdat-2*truncparam,1);
   
   % Define the filtered series (Note 2K observations are lost)
   for ii=truncparam+1:lvdat-truncparam
      filterseries1(ii-truncparam)=a'*series1(ii-truncparam:ii+truncparam);
      filterseries2(ii-truncparam)=a'*series2(ii-truncparam:ii+truncparam);
   end
   
   % Define the correlation for each periodicity
   tempcorr=corrcoef(filterseries1,filterseries2);
   corr(nn)=tempcorr(numvars,1);
end