% problem set 3 for asset pricing, longer vars 

clear all
close all

x = load('..\data\vwr.txt');

vw_dates = x(:,1); 
R  = 1+x(:,2); 
Rx = 1+x(:,3);
DP = R./Rx-1; % (D_T+P_t)/(P_t-1) * P_t-1/P_t
DD = [NaN; DP(2:end)./DP(1:end-1).*Rx(2:end)];  % D_t/P_t / (D_t-1/P_t-1) *P_t/P_t-1

T = size(R,1); 
r = log(R);
dp = log(DP); 
dd = log(DD) ; % can't take log of -99. make sure not to use first dd

disp('rho = mean D/P / (1+mean ( D/P))');
rho = mean(1./DP)/(1+mean(1./DP)); 
disp(rho)

ddi = [NaN; r(2:end)+rho*dp(2:end)-dp(1:end-1)];


% ----------------------------
% return and dp var. Include 1 more lag of dp than of other variables.
% you can also do a return and dd var, or both if you handle the identity
% right. 
% ----------------------------

lags = 0; % always dp-1, then say r-1, dp-1 and dp-2 = lag "1"
T = size(x,1);
rhv = [ones(T-1-lags,1) dp(1+lags:T-1)];
[b,seb,R2,R2vadj,v,F] = olsgmm([dp(1+lags+1:T) r(1+lags+1:T) dd(1+lags+1:T) ddi(1+lags+1:T)],rhv,0,0);
fit = rhv*b;

fprintf('var coeffs, \n')
fprintf('    '); fprintf( ' %6.6s ','dp_t'); for i = 1:lags; fprintf( ' %5.5s%1.0f ','dp_t-',i);fprintf( ' %5.5s%1.0f ','dd_t-',i); end; fprintf( ' %6.6s ','R2'); fprintf('\n');
fprintf(' dp '); fprintf( ' %6.3f ',b(2:end,1)); fprintf( '%6.3f ',R2(1)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,1)./seb(2:end,1)); fprintf('\n');
fprintf('  r '); fprintf( ' %6.3f ',b(2:end,2)); fprintf( ' %6.3f ',R2(2)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,2)./seb(2:end,2)); fprintf('\n');
fprintf(' dd '); fprintf( ' %6.3f ',b(2:end,3)); fprintf( ' %6.3f ',R2(3)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,3)./seb(2:end,3)); fprintf('\n');
fprintf(' ddi'); fprintf( ' %6.3f ',b(2:end,4)); fprintf( ' %6.3f ',R2(3)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,4)./seb(2:end,4)); fprintf('\n');

disp(' '); 

title1 = ['Dp shock, no \Deltad, full sample'];
title2 = ['Dd shock, no dp, full sample'];
var_and_plot; % script to run VAR and plot
print -depsc2 ps3_full_nolag.eps

% ----------------------------
% Lag 1 in full sample
% ----------------------------

% DD on rhs 
lags = 1; % always dp-1, then say r-1, dp-1 and dp-2 = lag "1"
T = size(x,1);
rhv = [ones(T-1-lags,1) dp(1+lags:T-1)];
for i = 1:lags
    rhv = [rhv dp(1+lags-i:T-1-i) ddi(1+lags-i+1:T-i)];
end;
[b,seb,R2,R2vadj,v,F] = olsgmm([dp(1+lags+1:T) r(1+lags+1:T) ddi(1+lags+1:T)],rhv,0,0);
fit = rhv*b;


fprintf('var coeffs, dp-1, dp-2, dd-1, dp-3, dd-2, etc. Last col is R2 \n')
fprintf('    '); fprintf( ' %6.6s ','dp_t'); for i = 1:lags; fprintf( ' %5.5s%1.0f ','dp_t-',i-1);fprintf( ' %5.5s%1.0f ','dd_t-',i-1); end; fprintf( ' %6.6s ','R2'); fprintf('\n');
fprintf(' dp '); fprintf( ' %6.3f ',b(2:end,1)); fprintf( '%6.3f ',R2(1)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,1)./seb(2:end,1)); fprintf('\n');
fprintf('  r '); fprintf( ' %6.3f ',b(2:end,2)); fprintf( ' %6.3f ',R2(2)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,2)./seb(2:end,2)); fprintf('\n');
fprintf('ddi '); fprintf( ' %6.3f ',b(2:end,3)); fprintf( ' %6.3f ',R2(3)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,3)./seb(2:end,3)); fprintf('\n');

disp(' '); 

% DD,R on rhs 
lags = 1; % always dp-1, then say r-1, dp-1 and dp-2 = lag "1"
T = size(x,1);
rhv = [ones(T-1-lags,1) dp(1+lags:T-1) r(1+lags:T-1) ddi(1+lags:T-1)];
[b,seb,R2,R2vadj,v,F] = olsgmm([dp(1+lags+1:T) r(1+lags+1:T) ddi(1+lags+1:T)],rhv,0,0);

fprintf('    '); fprintf( ' %6.6s ','dp_t'); fprintf( ' %6.6s ','r_t');fprintf( ' %6.6s ','dd_t'); fprintf( ' %6.6s ','R2'); fprintf('\n');
fprintf(' dp '); fprintf( ' %6.3f ',b(2:end,1)); fprintf( '%6.3f ',R2(1)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,1)./seb(2:end,1)); fprintf('\n');
fprintf('  r '); fprintf( ' %6.3f ',b(2:end,2)); fprintf( ' %6.3f ',R2(2)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,2)./seb(2:end,2)); fprintf('\n');
fprintf('ddi '); fprintf( ' %6.3f ',b(2:end,3)); fprintf( ' %6.3f ',R2(3)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,3)./seb(2:end,3)); fprintf('\n');

% now R on rhs
T = size(x,1);
rhv = [ones(T-1-lags,1) dp(1+lags:T-1)];
for i = 1:lags
    rhv = [rhv dp(1+lags-i:T-1-i) r(1+lags-i+1:T-i)];
end;
[b,seb,R2,R2vadj,v,F] = olsgmm([dp(1+lags+1:T) r(1+lags+1:T) ddi(1+lags+1:T)],rhv,0,0);
fit = rhv*b;

fprintf('var coeffs, dp-1, dp-2, r-1, dp-3, r-2, etc. Last col is R2 \n')
fprintf('    '); fprintf( ' %6.6s ','dp_t'); for i = 1:lags; fprintf( ' %5.5s%1.0f ','dp_t-',i-1);fprintf( ' %5.5s%1.0f ','r_t-',i-1); end; fprintf( ' %6.6s ','R2'); fprintf('\n');
fprintf(' dp '); fprintf( ' %6.3f ',b(2:end,1)); fprintf( '%6.3f ',R2(1)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,1)./seb(2:end,1)); fprintf('\n');
fprintf('  r '); fprintf( ' %6.3f ',b(2:end,2)); fprintf( ' %6.3f ',R2(2)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,2)./seb(2:end,2)); fprintf('\n');
fprintf('ddi '); fprintf( ' %6.3f ',b(2:end,3)); fprintf( ' %6.3f ',R2(3)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,3)./seb(2:end,3)); fprintf('\n');


title1 = ['Dp shock, no \Deltad, lags = ' num2str(lags, '%3.0f')];
title2 = ['Dd shock, no dp, lags = ' num2str(lags, '%3.0f')];
var_and_plot; % script to run VAR and plot

print -depsc2 ps3_one_lag_long.eps

% ----------------------
% many lags 
% ---------------------- 

lags = 2; 
rhv = [ones(T-1-lags,1) dp(1+lags:T-1)];
for i = 1:lags
    rhv = [rhv dp(1+lags-i:T-1-i) r(1+lags-i+1:T-i)];
end;
[b,seb,R2,R2vadj,v,F] = olsgmm([dp(1+lags+1:T) r(1+lags+1:T) ddi(1+lags+1:T)],rhv,0,0);

fprintf('var coeffs, dp-1, dp-2, r-1, dp-3, r-2, etc. Last col is R2 \n')
fprintf('    '); fprintf( ' %6.6s ','dp_t'); for i = 1:lags; fprintf( ' %5.5s%1.0f ','dp_t-',i-1);fprintf( ' %5.5s%1.0f ','r_t-',i-1); end; fprintf( ' %6.6s ','R2'); fprintf('\n');
fprintf(' dp '); fprintf( ' %6.3f ',b(2:end,1)); fprintf( '%6.3f ',R2(1)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,1)./seb(2:end,1)); fprintf('\n');
fprintf('  r '); fprintf( ' %6.3f ',b(2:end,2)); fprintf( ' %6.3f ',R2(2)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,2)./seb(2:end,2)); fprintf('\n');
fprintf('ddi '); fprintf( ' %6.3f ',b(2:end,3)); fprintf( ' %6.3f ',R2(3)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,3)./seb(2:end,3)); fprintf('\n');

title1 = ['Dp shock, no \Deltad, lags = ' num2str(lags, '%3.0f')];
title2 = ['Dd shock, no dp, lags = ' num2str(lags, '%3.0f')];
var_and_plot; % script to run VAR and plot
print -depsc2 ps3_resp2_lags2.eps;


lags = 3; 
rhv = [ones(T-1-lags,1) dp(1+lags:T-1)];
for i = 1:lags
    rhv = [rhv dp(1+lags-i:T-1-i) r(1+lags-i+1:T-i)];
end;
[b,seb,R2,R2vadj,v,F] = olsgmm([dp(1+lags+1:T) r(1+lags+1:T) ddi(1+lags+1:T)],rhv,0,0);
fit = rhv*b;

fprintf('var coeffs, dp-1, dp-2, r-1, dp-3, r-2, etc. Last col is R2 \n')
fprintf('    '); fprintf( ' %6.6s ','dp_t'); for i = 1:lags; fprintf( ' %5.5s%1.0f ','dp_t-',i-1);fprintf( ' %5.5s%1.0f ','r_t-',i-1); end; fprintf( ' %6.6s ','R2'); fprintf('\n');
fprintf(' dp '); fprintf( ' %6.3f ',b(2:end,1)); fprintf( '%6.3f ',R2(1)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,1)./seb(2:end,1)); fprintf('\n');
fprintf('  r '); fprintf( ' %6.3f ',b(2:end,2)); fprintf( ' %6.3f ',R2(2)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,2)./seb(2:end,2)); fprintf('\n');
fprintf('ddi '); fprintf( ' %6.3f ',b(2:end,3)); fprintf( ' %6.3f ',R2(3)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,3)./seb(2:end,3)); fprintf('\n');

title1 = ['Dp shock, no \Deltad, lags = ' num2str(lags, '%3.0f')];
title2 = ['Dd shock, no dp, lags = ' num2str(lags, '%3.0f')];
var_and_plot; % script to run VAR and plot
print -depsc2 ps3_resp2_lags3.eps;

lags = 4; 
rhv = [ones(T-1-lags,1) dp(1+lags:T-1)];
for i = 1:lags
    rhv = [rhv dp(1+lags-i:T-1-i) r(1+lags-i+1:T-i)];
end;
[b,seb,R2,R2vadj,v,F] = olsgmm([dp(1+lags+1:T) r(1+lags+1:T) ddi(1+lags+1:T)],rhv,0,0);

fprintf('var coeffs, dp-1, dp-2, r-1, dp-3, r-2, etc. Last col is R2 \n')
fprintf('    '); fprintf( ' %6.6s ','dp_t'); for i = 1:lags; fprintf( ' %5.5s%1.0f ','dp_t-',i-1);fprintf( ' %5.5s%1.0f ','r_t-',i-1); end; fprintf( ' %6.6s ','R2'); fprintf('\n');
fprintf(' dp '); fprintf( ' %6.3f ',b(2:end,1)); fprintf( '%6.3f ',R2(1)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,1)./seb(2:end,1)); fprintf('\n');
fprintf('  r '); fprintf( ' %6.3f ',b(2:end,2)); fprintf( ' %6.3f ',R2(2)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,2)./seb(2:end,2)); fprintf('\n');
fprintf('ddi '); fprintf( ' %6.3f ',b(2:end,3)); fprintf( ' %6.3f ',R2(3)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,3)./seb(2:end,3)); fprintf('\n');

title1 = ['Dp shock, no \Deltad, lags = ' num2str(lags, '%3.0f')];
title2 = ['Dd shock, no dp, lags = ' num2str(lags, '%3.0f')];
var_and_plot; % script to run VAR and plot
print -depsc2 ps3_resp2_lags4.eps;

% ----------------------
% Postwar
% ---------------------

st = 27; 
'sample' 
[vw_dates(st) vw_dates(end)]; 
dd = dd(st:end);
dp = dp(st:end);
r = r(st:end); 
ddi = ddi(st:end); 
T = size(dd,1);

% ------------------------
% zero lags for reference
% ------------------------

lags = 0; % always dp-1, then say r-1, dp-1 and dp-2 = lag "1"
rhv = [ones(T-1-lags,1) dp(1+lags:T-1)];
[b,seb,R2,R2vadj,v,F] = olsgmm([dp(1+lags+1:T) r(1+lags+1:T) dd(1+lags+1:T) ddi(1+lags+1:T)],rhv,0,0);
fit = rhv*b;

fprintf('var coeffs, \n')
fprintf('    '); fprintf( ' %6.6s ','dp_t'); for i = 1:lags; fprintf( ' %5.5s%1.0f ','dp_t-',i);fprintf( ' %5.5s%1.0f ','dd_t-',i); end; fprintf( ' %6.6s ','R2'); fprintf('\n');
fprintf(' dp '); fprintf( ' %6.3f ',b(2:end,1)); fprintf( '%6.3f ',R2(1)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,1)./seb(2:end,1)); fprintf('\n');
fprintf('  r '); fprintf( ' %6.3f ',b(2:end,2)); fprintf( ' %6.3f ',R2(2)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,2)./seb(2:end,2)); fprintf('\n');
fprintf(' dd '); fprintf( ' %6.3f ',b(2:end,3)); fprintf( ' %6.3f ',R2(3)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,3)./seb(2:end,3)); fprintf('\n');
fprintf(' ddi'); fprintf( ' %6.3f ',b(2:end,4)); fprintf( ' %6.3f ',R2(3)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,4)./seb(2:end,4)); fprintf('\n');

disp(' '); 

title1 = ['Dp shock, no \Deltad, postwar sample'];
title2 = ['Dd shock, no dp, postwar sample'];
var_and_plot; % script to run VAR and plot
print -depsc2 ps3_postwar_nolag.eps



% ------------------------
% 1 lags 
% ------------------------

%dd - dp
lags = 1; % always dp-1, then say r-1, dp-1 and dp-2 = lag "1"
rhv = [ones(T-1-lags,1) dp(1+lags:T-1)];
for i = 1:lags
    rhv = [rhv dp(1+lags-i:T-1-i) ddi(1+lags-i+1:T-i)];
end;
[b,seb,R2,R2vadj,v,F] = olsgmm([dp(1+lags+1:T) r(1+lags+1:T) ddi(1+lags+1:T)],rhv,0,0);

fprintf('var coeffs, dp-1, dp-2, r-1, dp-3, r-2, etc. Last col is R2 \n')
fprintf('    '); fprintf( ' %6.6s ','dp_t'); for i = 1:lags; fprintf( ' %5.5s%1.0f ','dp_t-',i-1);fprintf( ' %5.5s%1.0f ','dd_t-',i-1); end; fprintf( ' %6.6s ','R2'); fprintf('\n');
fprintf(' dp '); fprintf( ' %6.3f ',b(2:end,1)); fprintf( '%6.3f ',R2(1)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,1)./seb(2:end,1)); fprintf('\n');
fprintf('  r '); fprintf( ' %6.3f ',b(2:end,2)); fprintf( ' %6.3f ',R2(2)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,2)./seb(2:end,2)); fprintf('\n');
fprintf('ddi '); fprintf( ' %6.3f ',b(2:end,3)); fprintf( ' %6.3f ',R2(3)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,3)./seb(2:end,3)); fprintf('\n');

%dd - r
lags = 1; % always dp-1, then say r-1, dp-1 and dp-2 = lag "1"
rhv = [ones(T-1-lags,1) dp(1+lags:T-1)];
for i = 1:lags
    rhv = [rhv r(1+lags-i+1:T-i) ddi(1+lags-i+1:T-i)];
end;
[b,seb,R2,R2vadj,v,F] = olsgmm([dp(1+lags+1:T) r(1+lags+1:T) ddi(1+lags+1:T)],rhv,0,0);

fprintf('    '); fprintf( ' %6.6s ','dp_t'); for i = 1:lags; fprintf( ' %5.5s%1.0f ','r_t-',i-1);fprintf( ' %5.5s%1.0f ','dd_t-',i-1); end; fprintf( ' %6.6s ','R2'); fprintf('\n');
fprintf(' dp '); fprintf( ' %6.3f ',b(2:end,1)); fprintf( '%6.3f ',R2(1)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,1)./seb(2:end,1)); fprintf('\n');
fprintf('  r '); fprintf( ' %6.3f ',b(2:end,2)); fprintf( ' %6.3f ',R2(2)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,2)./seb(2:end,2)); fprintf('\n');
fprintf('ddi '); fprintf( ' %6.3f ',b(2:end,3)); fprintf( ' %6.3f ',R2(3)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,3)./seb(2:end,3)); fprintf('\n');

%dp-r
lags = 1; % always dp-1, then say r-1, dp-1 and dp-2 = lag "1"
rhv = [ones(T-1-lags,1) dp(1+lags:T-1)];
for i = 1:lags
    rhv = [rhv dp(1+lags-i:T-1-i) r(1+lags-i+1:T-i)];
end;
[b,seb,R2,R2vadj,v,F] = olsgmm([dp(1+lags+1:T) r(1+lags+1:T) ddi(1+lags+1:T)],rhv,0,0);

fprintf('var coeffs, dp-1, dp-2, r-1, dp-3, r-2, etc. Last col is R2 \n')
fprintf('    '); fprintf( ' %6.6s ','dp_t'); for i = 1:lags; fprintf( ' %5.5s%1.0f ','dp_t-',i-1);fprintf( ' %5.5s%1.0f ','r_t-',i-1); end; fprintf( ' %6.6s ','R2'); fprintf('\n');
fprintf(' dp '); fprintf( ' %6.3f ',b(2:end,1)); fprintf( '%6.3f ',R2(1)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,1)./seb(2:end,1)); fprintf('\n');
fprintf('  r '); fprintf( ' %6.3f ',b(2:end,2)); fprintf( ' %6.3f ',R2(2)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,2)./seb(2:end,2)); fprintf('\n');
fprintf('ddi '); fprintf( ' %6.3f ',b(2:end,3)); fprintf( ' %6.3f ',R2(3)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,3)./seb(2:end,3)); fprintf('\n');

title1 = ['Dp shock, no \Deltad, 1 lag postwar'];
title2 = ['Dd shock, no dp, 1 lag postwar'];
var_and_plot; % script to run VAR and plot
print -depsc2 ps3_postwar_1_lag.eps

% ---------------------
% produce new shocks and responses
% ----------------------

% now define new shocks that are linear combinations of the old ones (same
% as new ir that are linear combinations of old ones so that the response
% of sum rho^jr_t+j to shock 1 is zero and the response of sum rho^j dd_t+j
% to the other shock is zero

% if the e_dp (vector) shock produces sumr_dp and the e_dd shock(vector)
% produces sumr_dd then (1-a)*edp + a*edd produces 
% a*sumr_dp->r + b*sumr_dd->r = 1
% a*sumr_dp->d + b*sumr_dd->d.= 0

%[ sumr_dp->r + sumr_dd->r ] [a]    [1] 
%[ sumr_dp->d + sumr_dd->d.] [b] =  [0]

% needs det > 0 which isn't the case for the 0 lag model. 

% define one shock with sum rho ^1 r_1 + rho^2 r_2 ... while the other with
% sum rho^0 d0 + rho ^1 d1 etc. These are the terms of the return
% decomposition

% xv_dp = [dp dp-1 r cumr_dp ddt_dp dt_dp chpt_dp pt_dp]; 

rj = (rho.^(0:1:size(xv,1)-2)'*ones(1,4)).*[xv_dp(2:end,[3 5]) xv_dd(2:end,[3 5])];
sumr_dp = sum(rj(2:end,1));
sumr = sum(rj(2:end,3));
sumd_dp = sum(rj(1:end,2)); % returns from 2 on, all dividends. 
sumd = sum(rj(1:end,4));

[sumr_dp sumr sumd_dp sumd]

mat = [sumr_dp sumr ; sumd_dp sumd];
if abs(det(mat))> 0; 
    abr = inv(mat)*[1 ; 0];
    abd = inv(mat)*[0 ; 1];
    
    xv_r = xv_dp*abr(1)+ xv_dd*abr(2); 
    xv_d = xv_dp*abd(1)+ xv_dd*abd(2);
    
    figure
    subplot(1,2,1);
    xaxis = (1:pmax)';
    plot(xaxis, xv_r(1:pmax,1), '-g.','linewidth',2)
    hold on
    plot(xaxis, xv_r(1:pmax,3), '-b.','linewidth',2);
    plot(xaxis, xv_r(1:pmax,4), '-ob','linewidth',2);
    plot(xaxis, xv_r(1:pmax,8), '-k.','linewidth',2);
    plot(xaxis, xv_r(1:pmax,5),'-r.','linewidth',2);
    plot(xaxis, xv_r(1:pmax,6),'-vr','linewidth',2);
    plot(xaxis,0*xaxis,'-k')
    legend('dp','r','cumr','p','\Deltad','d');
    
    sumr = sum(rho.^(1:ir_T-2)'.*xv_r(3:end,3)); % sum up only future terms, not current
    sumd = sum(rho.^(0:ir_T-2)'.*xv_r(2:end,5));
    text(15,-0.4,['\Sigma_{j=2}\rho^{j-1}r_{j}=' num2str(sumr,'%4.2f')],'color','b')
    text(15,-0.6,['\Sigma_{j=1}\rho^{j-1}\Deltad_{j}=' num2str(sumd,'%4.2f')],'color','r')
    
    %axis([0 10 -0.2 1.2]);
    title(['long run r shock lags = ' num2str(lags, '%3.0f')]);
    
    
    subplot(1,2,2);
    xaxis = (1:pmax)';
    plot(xaxis, xv_d(1:pmax,1), '-g.','linewidth',2)
    hold on
    plot(xaxis, xv_d(1:pmax,3), '-b.','linewidth',2);
    plot(xaxis, xv_d(1:pmax,4), '-ob','linewidth',2);
    plot(xaxis, xv_d(1:pmax,8), '-k.','linewidth',2);
    plot(xaxis, xv_d(1:pmax,5),'-r.','linewidth',2);
    plot(xaxis, xv_d(1:pmax,6),'-vr','linewidth',2);
    plot(xaxis,0*xaxis,'-k')
    legend('dp','r','cumr','p','\Deltad','d');
    
    sumr = sum(rho.^(1:ir_T-2)'.*xv_d(3:end,3)); % sum up only future terms, not current
    sumd = sum(rho.^(0:ir_T-2)'.*xv_d(2:end,5)); % includes current term
    text(15,0.5,['\Sigma_{j=2}\rho^{j-1}r_{j}=' num2str(sumr,'%4.2f')],'color','b')
    text(15,0.3,['\Sigma_{j=1}\rho^{j}\Deltad_{j}=' num2str(sumd,'%4.2f')],'color','r')
    
    %axis([0 10 -0.2 1.2]);
    title(['long run d shock lags = ' num2str(lags, '%3.0f')]);
    
    set(gcf, 'PaperPositionMode', 'manual');
    set(gcf, 'PaperUnits', 'inches');
    set(gcf, 'PaperPosition', [0.25 2.5 10.0 4.0]);
    
    print -depsc2 ps3_newshock_lag1.eps

    disp(' ');
    disp('time 0 responses -- dp r dd  '); 
    disp( xv_r(2,[1 3 5] ));
    disp( xv_d(2,[1 3 5] ));
        
end;

stop

% -------------
% Two lags postwar
% ---------------
lags = 2; % always dp-1, then say r-1, dp-1 and dp-2 = lag "1"
rhv = [ones(T-1-lags,1) dp(1+lags:T-1)];
for i = 1:lags
    rhv = [rhv dp(1+lags-i:T-1-i) r(1+lags-i+1:T-i)];
end;
[b,seb,R2,R2vadj,v,F] = olsgmm([dp(1+lags+1:T) r(1+lags+1:T) ddi(1+lags+1:T)],rhv,0,0);

fprintf('var coeffs, dp-1, dp-2, r-1, dp-3, r-2, etc. Last col is R2 \n')
fprintf('    '); fprintf( ' %6.6s ','dp_t'); for i = 1:lags; fprintf( ' %5.5s%1.0f ','dp_t-',i-1);fprintf( ' %5.5s%1.0f ','r_t-',i-1); end; fprintf( ' %6.6s ','R2'); fprintf('\n');
fprintf(' dp '); fprintf( ' %6.3f ',b(2:end,1)); fprintf( '%6.3f ',R2(1)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,1)./seb(2:end,1)); fprintf('\n');
fprintf('  r '); fprintf( ' %6.3f ',b(2:end,2)); fprintf( ' %6.3f ',R2(2)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,2)./seb(2:end,2)); fprintf('\n');
fprintf('ddi '); fprintf( ' %6.3f ',b(2:end,3)); fprintf( ' %6.3f ',R2(3)); fprintf('\n');
fprintf('  t '); fprintf( '(%6.3f)',b(2:end,3)./seb(2:end,3)); fprintf('\n');

title1 = ['Dp shock, no \Deltad, 2 lag postwar'];
title2 = ['Dd shock, no dp, 2 lag postwar'];
var_and_plot; % script to run VAR and plot
print -depsc2 ps3_postwar_2_lag.eps

