% ps3_time_series
% programs for problem set 3 time series component


clear all; 
close all; % good to start by clearing out stuff from previous runs

rows = inline('size(x,1)'); % defines shorter function to return size of matrices
cols = inline('size(x),2)'); 

addpath('../Data');  % modify according to your directory structure. This is where to find regress_jc on my computer
addpath('../Programs');

format short; % convenient output formats
format compact; 

 % *******************************
 % VAR with dividend growth and dp 
 % *******************************

load ps1_data.txt; 
cal = ps1_data(:,1); 
r = ps1_data(:,2); 
dp = ps1_data(:,3); 
dd = ps1_data(:,4); 
tbret = ps1_data(:,5); 

T = rows(r); 
maxdate = 1925+T; 

lnr = log(1+r);
lndp = log(dp); 
lndd = [-99; log(1+dd(2:end))] ; % can't take log of -99. make sure not to use first dd

rho = 0.96; 

% regressions on dp

b = [ones(T-1,1) lndp(1:T-1)]\[lnr(2:T) lndd(2:T) lndp(2:T)]; 
% this is a compact way to run three  regressions at once 
% b is now [ r constant   dd constant   dp constant ] 
%          [ r coeff      dd coeff      dp coeff    ] 

br = b(2,1);
bd = b(2,2);
bdp = b(2,3); 

disp('regression coefficients br bd and bdp'); 
disp([br bd bdp]); 

% keep errors for later
e = [lnr(2:T) lndd(2:T) lndp(2:T)]-[ones(T-1,1) lndp(1:T-1)]*b; 
er = e(:,1); 
ed = e(:,2); 
edp = e(:,3); 

% ------------------------------------------------------
% price and return decomposition
% -------------------------------------------------------

dlr = bd/(1-rho*bdp); 
rlr = br/(1-rho*bdp); 

lndp1 = lndp-mean(lndp); 
lndd1 = lndd-mean(lndd(2:end));
lnr1 = lnr-mean(lnr(2:end));

figure; 
plot((1926:maxdate)',[-lndp1 dlr*lndp1 -rlr*lndp1],'linewidth',2); 
legend('p-d','dividend component','return component',2);
print -depsc2 ps2_dp_decomp.eps; 

% ex post rational

T = size(lndp1,1); 
rterm = NaN(T,1);
dterm = NaN(T,1);
pdterm = NaN(T,1);
for t = 1:T-1; 
    rterm(t) = sum((rho.^(0:1:T-t-1)').*lnr1(t+1:T));
    dterm(t) = sum((rho.^(0:1:T-t-1)').*lndd1(t+1:T));
    pdterm(t) = rho^(T-t-1)*(-lndp1(T)); 
end;
rterm(T)=0;
dterm(T)=0;
pdterm(T)=-lndp1(T); 

figure;
plot((1926:maxdate)',[-lndp1 dterm-rterm+pdterm],'linewidth',2); 
legend('p-d','identity',2);
print -depsc2 check_id.eps

figure;
%plot((1926:maxdate)',[dterm-rterm+pdterm dterm pdterm],'linewidth',2); 
plot((1926:maxdate)',[-lndp1 dterm pdterm],'linewidth',2); 
hold on
plot((1926:maxdate)',0*(1926:maxdate)','-k'); 
legend('p-d','dividend component','terminal pd',2);
print -depsc2 ps2_dp_decomp_expost_d.eps; 

figure;
%plot((1926:maxdate)',[dterm-rterm+pdterm -rterm pdterm]); 
plot((1926:maxdate)',[-lndp1 -rterm pdterm],'linewidth',2); 
hold on
plot((1926:maxdate)',0*(1926:maxdate)','-k'); 
legend('p-d','return component','terminal pd',2);
print -depsc2 ps2_dp_decomp_expost_r.eps; 


% ------------------
% Now, return decomposition 
% -------------------------

lnr1 = lnr - mean(lnr); 
figure; 
plot((1927:maxdate)',[lnr1(2:end) br*lndp1(1:T-1)],'linewidth',2); 
legend('r_t','E_{t-1}(r_t)'); 
axis([1927 maxdate -inf inf]); 
print -depsc2 ps2_r_decomp1.eps; 

dcoef2 = bd*rho/(1-rho*bdp); 
rcoef2 = br*rho/(1-rho*bdp); 

figure; 
plot((1927:maxdate)',[ ed+0.3 dcoef2*edp -rcoef2*edp],'linewidth',2); 
legend('\Deltad_t term+0.3','future \Deltad','future r',4); 
axis([1927 maxdate -inf inf]);
print -depsc2 ps2_r_decomp2.eps; 

figure; 
plot((1927:maxdate)',[er ed+dcoef2*edp-rcoef2*edp],'linewidth',2);
print -depsc2 ps2_r_decomp_check_id.eps;
% check that terms add up

% ---------------------
% variance decomposition
% ---------------------

disp('100 x var of r, cov(r,d), cov(r,futured) cov(r,futurer)');
crd=cov(er,ed);
crd = crd(1,2);
crdp = cov(er,edp);
crdp = crdp(1,2);
disp(100*[var(er) crd rho/(1-rho*bdp)*bd*crdp -rho/(1-rho*bdp)*br*crdp]);
disp('as percent var (r)'); 
disp(100*[var(er) crd rho/(1-rho*bdp)*bd*crdp -rho/(1-rho*bdp)*br*crdp]./var(er));

