% ps2
% programs for problem set 2


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

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; 

% -----------------------------------------
% switches to control how the program runs
% -----------------------------------------

printdata = 1; 
display = 0; 

% *******************************
% 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 = size(r,1); 
maxdate = 1925+T; 
disp('data sample')
disp([cal(1,:) cal(end,:)]);

lnr = log(1+r);
lndp = log(dp); 
lndd = log(1+dd) ; % Make sure not to use first dd

rho = 0.96; 

% --------------------------
% class question. can we recover rho and the identity by regression? 
% ------------------------------
if 0;
    
    T = size(lnr,1);
    %lnr_t+1 = k - rho*lndp_t+1 + lndd_t+1 + lndp_t
    %lnr_t+1 - lndd_t+1 - lndp_t = k - rho*lndp_t+1
    
    y = lnr(2:T)-lndd(2:T)-lndp(1:T-1);
    x = [ones(T-1,1) lndp(2:T)];
    b = x\y;
    k = b(1);
    rho = -b(2);
    R2 = var(x*b)/var(y);
    R2R = 1-var(lnr(2:T)-(k-rho*lndp(2:T)+lndd(2:T)+lndp(1:T-1)))/var(lnr(2:T));
    disp('k and rho from regression; R2 R2(r)')
    disp([k rho R2 R2R]);
    
    %but why not lnr_t+1 = k - rho*lndp_t+1 + lndd_t+1 + lndp_t
    %lndp_t+1 = 1/rho k - 1/rho(lnr_t+1 - lndd_t+1 - lndp_t)
    
    x = [ones(T-1,1) lnr(2:T)-lndd(2:T)-lndp(1:T-1)];
    y = lndp(2:T);
    b = x\y;
    R2 = var(x*b)/var(y);
    rho = -1/b(2);
    k = rho*b(1);
    R2R = 1-var(lnr(2:T)-(k-rho*lndp(2:T)+lndd(2:T)+lndp(1:T-1)))/var(lnr(2:T));
    disp('k and rho from inverse regression;R2, R2r')
    disp([k rho R2 R2R]);
    
    PD = exp(mean(-lndp));
    DP = exp(mean(-lndp));
    rho = PD/(1+PD);
    disp('PD DP and rho, conventional');
    disp([PD DP rho]);
end;

%-------------------------- 
% create dividend growth by identity
% -------------------------

lndd_id = [NaN; lnr(2:T)+rho*lndp(2:T)-lndp(1:T-1)]; 

figure
plot((1926:maxdate)', [lndd-mean(lndd(2:end)) lndd_id-mean(lndd_id(2:end))])
title('Real and approximate identity dividend growth.')
print -depsc2 approxid.eps;

% regressions on dp

b = [ones(T-1,1) lndp(1:T-1)]\[lnr(2:T) lndd_id(2:T) lndp(2:T) ]; 
% this is a compact way to run three regressions at once 
% b is now [ dd constant   dp constant ] 
%          [ dd coeff      dp coeff    ] 
% Note you really should run all variables  on  dd,r too, but the coefficients come out near zero so I simplfy by ignoring that
br  = b(2,1);
bdd = b(2,2);
bdp = b(2,3); 


disp('regression coefficients of log r log dd and log dp on lagged log dp'); 
disp([ br bdd bdp]); 
disp('identity br = 1-rho*bdp +bd?')
disp([br  1-rho*bdp+bdd]); 

% keep errors for later
e = [lnr(2:T) lndd_id(2:T) lndp(2:T)]-[ones(T-1,1) lndp(1:T-1)]*b; 
er = e(:,1);
edd = e(:,2); 
edp = e(:,3); 
disp('first five r residuals, direct and implied') 
disp([er(1:5) -rho*edp(1:5)+edd(1:5)])

% using real dividend growth for comparison

b2 = [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 [ dd constant   dp constant ] 
%          [ dd coeff      dp coeff    ] 
% Note you really should run all variables  on  dd,r too, but the coefficients come out near zero so I simplfy by ignoring that
b2r  = b2(2,1);
b2dd = b2(2,2);
b2dp = b2(2,3); 


disp('regression coefficients of log r log dd and log dp on lagged log dp using actual dividend growth'); 
disp([ b2r b2dd b2dp]); 
disp('identity br = 1-rho*bdp +bd?')
disp([b2r  1-rho*b2dp+b2dd]); 

% now do the simulations
 
dpt = zeros(25,1); % set them up and leave the first element to be zero. 
ddt = zeros(25,1); 
rt = zeros(25,1); 
dpt(2) = 1; % first shock
ddt(2) = 0;
rt(2) = -rho; 

for t=3:25; % feed through the system
    dpt(t) = bdp*dpt(t-1);
    ddt(t) = bdd*dpt(t-1); 
    rt(t) = br*dpt(t-1);
end; 
dt = cumsum(ddt); % level of dividends is sum of changes

% use identities to get other variables
chpt = [0 ; (-dpt(2:end) + ddt(2:end) + dpt(1:end-1))];
pt = cumsum(chpt); 
cumrt = cumsum(rt);

figure;
xaxis = (1:25)';
plot(xaxis, dpt, '-', xaxis, cumrt, '-o',xaxis, pt, '-v',xaxis, dt,'-^'); 
legend('dp','\Sigmar','p','d'); 
hold on
plot(xaxis,0*xaxis,'-k');
title('Cumulative responses to a dp shock'); 
print -depsc2 ps2_erresp.eps; 

figure;
xaxis = (1:25)';
plot(xaxis, rt, '-o',xaxis, chpt, '-v',xaxis,ddt,'-^'); 
legend('r','\Deltap','\Deltad'); 
hold on
plot(xaxis,0*xaxis,'-k');
title('Raw responses to a dp shock'); 
print -depsc2 ps2_erresp_ch.eps; 

figure;
xaxis = (1:25)';
plot(xaxis, rt, '-o',xaxis, chpt, '-v',xaxis,ddt,'-^'); 
legend('r','\Deltap','\Deltad'); 
hold on
plot(xaxis,0*xaxis,'-k');
title('Raw responses to a dp shock'); 
axis([0 25 -0.05 0.15]);
print -depsc2 ps2_erresp_ch_big.eps; 

% second shock

rt = zeros(25,1); 
dpt = zeros(25,1); % set them up and leave the first element to be zero. 
ddt = zeros(25,1);
rt(2) = 1; 
dpt(2) = 0; % second shock
ddt(2) = 1;  

for t=3:25; % feed through the system
    dpt(t) = bdp*dpt(t-1);
    ddt(t) = bdd*dpt(t-1);
    rt(t) = br*dpt(t-1); 
end; 

dt = cumsum(ddt); % level of dividends is sum of changes

% use identities to get other variables
chpt = [0 ; (-dpt(2:end) + ddt(2:end) + dpt(1:end-1))];
pt = cumsum(chpt); 


figure; 
xaxis = (1:25)';
plot(xaxis, dpt, '-', xaxis, rt, '-o',xaxis, pt, '-v',xaxis, dt,'-^'); 
legend('dp','r','p','d'); 
hold on
plot(xaxis,0*xaxis,'-k');
axis([0 10 -0.2 1.2]); 
title('Response to a "dividend" shock'); 
print -depsc2 ps2_divresp.eps; 

% -----------------------
% what it ``should'' have looked like
% ----------------------------

br_old = 0;  
bdp_old = bdp;
bdd_old = br_old-1+rho*bdp_old;

disp('using coefficients to reflect 70s view: br bd bdp'); 
disp([br_old bdd_old bdp_old]); 

% now do the simulations -- same code with only coefficients changed
 
dpt = zeros(25,1); % set them up and leave the first element to be zero. 
ddt = zeros(25,1); 
rt = zeros(25,1); 
dpt(2) = 1; % first shock
ddt(2) = 0;
rt(2) = -rho; 

for t=3:25; % feed through the system
    dpt(t) = bdp_old*dpt(t-1);
    ddt(t) = bdd_old*dpt(t-1); 
    rt(t) = br_old*dpt(t-1);
end; 
dt = cumsum(ddt); % level of dividends is sum of changes

% use identities to get other variables
chpt = [0 ; (-dpt(2:end) + ddt(2:end) + dpt(1:end-1))];
pt = cumsum(chpt); 
cumrt = cumsum(rt);

figure;
xaxis = (1:25)';
plot(xaxis, dpt, '-', xaxis, cumrt, '-o',xaxis, pt, '-v',xaxis, dt,'-^'); 
legend('dp','\Sigmar','p','d'); 
hold on
plot(xaxis,0*xaxis,'-k');
axis([0 25 -1.6 1.1]);
title('Cumulative responses to a dp shock'); 
print -depsc2 ps2_erresp_old.eps; 

figure;
xaxis = (1:25)';
plot(xaxis, rt, '-o',xaxis, chpt, '-v',xaxis,ddt,'-^'); 
legend('r','\Deltap','\Deltad'); 
hold on
plot(xaxis,0*xaxis,'-k');
axis([0 25 -1.1 0.1]); 
title('Raw responses to a dp shock'); 
print -depsc2 ps2_erresp_ch_old.eps; 



% second shock

rt = zeros(25,1); 
dpt = zeros(25,1); % set them up and leave the first element to be zero. 
ddt = zeros(25,1);
rt(2) = 1; 
dpt(2) = 0; % second shock
ddt(2) = 1;  

for t=3:25; % feed through the system
    dpt(t) = bdp_old*dpt(t-1);
    ddt(t) = bdd_old*dpt(t-1);
    rt(t) = br_old*dpt(t-1); 
end; 

dt = cumsum(ddt); % level of dividends is sum of changes

% use identities to get other variables
chpt = [0 ; (-dpt(2:end) + ddt(2:end) + dpt(1:end-1))];
pt = cumsum(chpt); 


figure; 
xaxis = (1:25)';
plot(xaxis, dpt, '-', xaxis, rt, '-o',xaxis, pt, '-v',xaxis, dt,'-^'); 
legend('dp','r','p','d'); 
hold on
plot(xaxis,0*xaxis,'-k');
axis([0 10 -0.2 1.2]); 
title('Response to a "dividend" shock'); 
print -depsc2 ps2_divresp_old.eps; 


% -----------------------------
% simulation using Er, Ed models 
% Note: longer version in ps2_10.m with 
% dividend growth world too. 
% --------------------------------

randn('state',5)
phi = 0.94; 
rho = 0.96; 
sigEr = 0.05; % std. dev. of expected returns. 
sigr = 0.18; % std dev of actual returns  
sigd = 0.10; % std dev of dividend growth. 

sigepsx = (1-phi^2)^0.5*sigEr;     % std dev of shocks to expected returns x. 
                                   % We want to deliver the desired std. dev. expected returns x, 
                                   % and sigma^2(x) = (1/(1-phi^2))sigma^2(epsx)
                                  
sigepsr = (sigr^2 - sigEr^2)^0.5;  % sigma^2(r) = sigma^2(Er) + sigma^2(epsilon) so this tells you what sigma(epsilon) to assume
sigepsd = (sigd^2)^0.5; 

disp('assumed sigma(epsilon_x), sigma(epslion_r), sigma(epsilon_d) '); 
disp([sigepsx sigepsr sigepsd]);
% note I gave you these directly. As you can see, I reverse engineered them
% to get the input values I wanted. 

T = 10000; 
T1 = 80; 
x = zeros(T,1); 
y = zeros(T,1); 
r = zeros(T,1);
dd = zeros(T,1); 
for t = 2:T; 
    x(t) = phi*x(t-1)+sigepsx*randn(1);
    r(t) = x(t-1)+sigepsr*randn(1);
    y(t) = -x(t); 
    dd(t) = sigepsd*randn(1); 
end; 

dp = 1./(1-rho*phi)*x ; 

figure; 
subplot(2,1,1)
plot((1:T1)', 100*[0; x(1:T1-1)], '-r','linewidth',2); 
axis([-inf inf -40 40]); 
hold on;
plot((1:T1)', 100*r(1:T1), '-k','linewidth',2); 
title('actual and expected return');  
subplot(2,1,2)
plot((1:T1)', 0*(1:T1)', '-r','linewidth',2); 
hold on
plot((1:T1)', 100*dd(1:T1), '-k','linewidth',2);
axis([-inf inf -20 20]); 
title('actual and expected dividend growth') 
print -depsc2 ps2_randEr.eps;


figure
plot((1:T1)',dp(1:T1), 'linewidth',2);
ylabel('Log D/P'); 
print -depsc2 ps2_dp_er.eps;

% regressions copied from first problem set
% changed to log r on log dp

disp(''); 
fprintf('Regression of log returns on lagged d-p in simulated data with %6.0f observations\n', T); 
fprintf('                  %10s %10s %10s %10s \n','b','se','t(b)','R2'); 

rhv = [ones(T-1,1) dp(1:T-1)]; 
lhv = r(2:end); 
[b,bse,bt,r2,F]  = regress_jc(lhv,rhv,display); 
fprintf('1 Yr.             %10.3f %10.3f %10.3f %10.3f \n', [b(2) bse(2) bt(2) r2]);

r5 = (r(1:T-4))+(r(2:T-3))+(r(3:T-2))+(r(4:T-1))+(r(5:T)); 
r5a = r5(2:5:size(r5,1)); 
dp5a = dp(1:5:(size(dp,1)-5)); 
rhv = [ones(size(dp5a,1),1) dp5a]; 
lhv = r5a; 
[b,bse,bt,r2,F]  = regress_jc(lhv,rhv,display); 
fprintf('5 Yr. no overlap  %10.3f %10.3f %10.3f %10.3f \n', [b(2) bse(2) bt(2) r2]);

% now dividends 
 
disp(''); 
fprintf('Regression of log dividend growth on lagged d-p in simulated data with %6.0f observations\n', T); 
fprintf('                  %10s %10s %10s %10s \n','b','se','t(b)','R2'); 

rhv = [ones(T-1,1) dp(1:T-1)]; 
lhv = dd(2:end); 
[b,bse,bt,r2,F]  = regress_jc(lhv,rhv,display); 
fprintf('1 Yr.             %10.3f %10.3f %10.3f %10.3f \n', [b(2) bse(2) bt(2) r2]);

dd5 = (dd(1:T-4))+(dd(2:T-3))+(dd(3:T-2))+(dd(4:T-1))+(dd(5:T)); 

dd5a = dd5(2:5:size(r5,1)); 
dp5a = dp(1:5:(size(dp,1)-5)); 
rhv = [ones(size(dp5a,1),1) dp5a]; 
lhv = dd5a; 
[b,bse,bt,r2,F]  = regress_jc(lhv,rhv,display); 
fprintf('5 Yr. no overlap  %10.3f %10.3f %10.3f %10.3f \n', [b(2) bse(2) bt(2) r2]);
