clc; format short g; clear all; close all;


% import data from Excel
X = importdata('econ7720_corr_data.xls')


% define variables
Y = X.data.Sheet1(:,1);                 % Y = Real GDP
P = X.data.Sheet1(:,2);                 % P = Real GDP Deflator    
M0 = X.data.Sheet2(:,1);                % MO = currency + reserves (monetary base)
M1 = X.data.Sheet2(:,2);                % M1 = currency + demand deposits
M2 = X.data.Sheet2(:,3);                % M2 = M1 + time deposits
R_FF = X.data.Sheet2(:,4);              % R_FF = Federal Funds Rate
R_3M = X.data.Sheet2(:,5);              % R_3M = 3-month Treasury Bill Rate
R_1Y = X.data.Sheet2(:,6);              % R_1Y = 1-year Treasury Bill Rate
R_10Y = X.data.Sheet2(:,7);             % R_10Y = 10-year Treasury Bill Rate


% transform monthly series into quarterly averages
Xq = mth2qtr([M0 M1 M2 R_FF R_3M R_1Y R_10Y],0);

M0 = Xq(:,1);
M1 = Xq(:,2);
M2 = Xq(:,3);
R_FF = Xq(:,4);
R_3M = Xq(:,5);
R_1Y = Xq(:,6);
R_10Y = Xq(:,7);


% log and detrend {Y,P,M0,M1,M2}
[smooth,deviation]=hpfilter([log(Y) log(P) log(M0) log(M1) log(M2)],1600);

y = deviation(:,1);                     % y = deviation of log Real GDP from trend
p = deviation(:,2);                     % p = deviation of log Real GDP Deflator from trend
m0 = deviation(:,3);                    % m0 = deviation of log M0 from trend
m1 = deviation(:,4);                    % m1 = deviation of log M1 from trend
m2 = deviation(:,5);                    % m2 = deviation of log M2 from trend

 
% compute sample autocorrelations
output1 = sample_moments([y m0 m1 m2],8)
output2 = sample_moments([y p R_FF R_1Y R_10Y],8)

corr_ym0 = [reshape(output1.autocorr(1,2,9:-1:1),1,9) reshape(output1.autocorr(2,1,2:1:9),1,8)];
corr_ym1 = [reshape(output1.autocorr(1,3,9:-1:1),1,9) reshape(output1.autocorr(3,1,2:1:9),1,8)];
corr_ym2 = [reshape(output1.autocorr(1,4,9:-1:1),1,9) reshape(output1.autocorr(4,1,2:1:9),1,8)];

corr_yp = [reshape(output2.autocorr(1,2,9:-1:1),1,9) reshape(output2.autocorr(2,1,2:1:9),1,8)];
corr_yff = [reshape(output2.autocorr(1,3,9:-1:1),1,9) reshape(output2.autocorr(3,1,2:1:9),1,8)];
corr_y1y = [reshape(output2.autocorr(1,4,9:-1:1),1,9) reshape(output2.autocorr(4,1,2:1:9),1,8)];
corr_y10y = [reshape(output2.autocorr(1,5,9:-1:1),1,9) reshape(output2.autocorr(5,1,2:1:9),1,8)];


% plot the sample autocorrelations
qtr = [-8:1:8];

figure,
plot(qtr,corr_ym0,'b',qtr,corr_ym1,'k',qtr,corr_ym2,'r'), grid on
title('Dynamic Correlations: corr(GDP_t,M_{t+j})')
xlabel('Quarters (j)')
legend('M0','M1','M2')

figure,
plot(qtr,corr_yp,'b',qtr,corr_yff,'g',qtr,corr_y1y,'r',qtr,corr_y10y,'k'), grid on
title('Dynamic Correlations: corr(Y_t,P_{t+j}) and corr(Y_t,R_{t+j})')
xlabel('Quarters (j)')
legend('P','FF','1Y','10Y') 


% reproduce time series plots from Walsh (page 16 and page 18)
time = [1967:.25:2004.75];

figure,
plot(time,100*y,'k-',time,100*m1,'k-.',time,100*m2,'k:')
title('Detrended Money and Real GDP')
legend('Y','M1','M2')
xlabel('year')
ylabel('percent')

figure,
plot(time,100*y,'k-',time,R_FF,'k-.',time,R_3M,'k:')
title('Interest Rates and Detrended Real GDP')
legend('Y','FF','3MTB')
xlabel('year')
ylabel('percent')