% problem set 2 for asset pricing, forecasting returns
% Simple VAR, univariate 

close all
clear 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); 

% ------------------------------------------------
% 0. Do impluse-response function in full data set
% -------------------------------------------------

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

disp('rh0 = mean D/P / (1+mean ( D/P))');
rho = mean(1./DP)/(1+mean(1./DP)); 
disp('rho from mean dp') 
disp(rho)
ron = 0.96 ; 
disp('using rho = 0.96');

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

figure; 
plot( floor(vw_dates(2:end)/10000) ,[dd(2:end)-mean(dd(2:end)) ddi(2:end)-mean(ddi(2:end))]); 
title('actual and implied dividend growth'); 
print -depsc2 ddplot

% -----------------------------
% un restricted VAR
% -----------------------------

T = size(r,1); 
[b,seb,R2,R2adj,v,F] = olsgmm([r(2:T) ddi(2:T) dd(2:T) dp(2:T) ],[ones(T-1,1) dp(1:T-1)],0,0);     

br = b(2,1);
bdi = b(2,2);
bd = b(2,3);
phi = b(2,4);
                                                     
fprintf('VAR coefficients %10.9s  %10.9s \n','dp_t-1','s.e'); 
fprintf('r                %10.3f (%10.3f)\n',br,seb(2,1));
fprintf('dd               %10.3f (%10.3f)\n',bd,seb(2,2));
fprintf('ddi              %10.3f (%10.3f)\n',bdi,seb(2,3));
fprintf('dp               %10.3f (%10.3f)\n',phi,seb(2,4));

fprintf('1-rho phi        %10.3f \n',1-rho*phi);

% error correlation

e = [r(2:T) ddi(2:T) dd(2:T) dp(2:T) ] - [ones(T-1,1) dp(1:T-1)]*b;
ecov = cov(e);
ecorr = corr(e); 
for i = 1:size(ecorr,1);
    ecorr(i,i) = ecov(i,i)^0.5;
end;

fprintf(' correlation / std dev matrix\n'); 
fprintf('         %10.9s %10.9s %10.9s %10.9s\n','r','ddi','dd','dp'); 
fprintf('r        %10.3f %10.3f %10.3f %10.3f \n',ecorr(1,:));
fprintf('ddi      %10.3f %10.3f %10.3f %10.3f \n',ecorr(2,:));
fprintf('dd       %10.3f %10.3f %10.3f %10.3f \n',ecorr(3,:));
fprintf('dp       %10.3f %10.3f %10.3f %10.3f \n',ecorr(4,:));


% var 

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 is epsilon x, shock to dp
rt(2) = -rho;

for t=3:25; % feed through the system
    dpt(t) = phi*dpt(t-1);
    ddt(t) = bdi*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

cumr = cumsum(rt);
chpt = [0 ; (-dpt(2:end) + ddt(2:end) + dpt(1:end-1))];
pt = cumsum(chpt); 

figure;
xaxis = (1:25)';
plot(xaxis, dpt, '-', xaxis, rt, '-r',xaxis,cumr,'-or',xaxis, pt, '-v',xaxis, ddt,'-k',xaxis, dt,'-^k'); 
hold on
plot(xaxis,0*xaxis,'-k')
legend('dp','r','sum r','p','\Deltad','d'); 
title('Response to a dp = "Er" shock'); 
print -depsc2 ps2_resp1a.eps; 


figure;
xaxis = (1:25)';
plot(xaxis, rt, '-r',xaxis,ddt,'-^k'); 
hold on
plot(xaxis,0*xaxis,'-k')
legend('r','\Deltad'); 
title('Response to a dp = "Er" shock'); 

axis([0 25 -0.15 0.15])
print -depsc2 ps2_resp1b.eps; 


% second shock

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

for t=3:25; % feed through the system
    dpt(t) = phi*dpt(t-1);
    ddt(t) = bd*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
% r_t+1 = d_t+1 - d_t - rho ( d_t+1 - p_t+1) + (d_t - p_t) 
cumr = cumsum(rt);
% p_t+1 - p_t = -(d_t+1 - p_t+1) + (d_t+1 - d_t ) + (d_t - p_t) 
chpt = [0 ; (-dpt(2:end) + ddt(2:end) + dpt(1:end-1))];
pt = cumsum(chpt); 

figure; 
xaxis = (1:25)';
plot(xaxis, dpt, '-', xaxis, rt, '-r',xaxis,cumr,'-or',xaxis, pt, '-v',xaxis, ddt,'-k',xaxis, dt,'-^k'); 
legend('dp','r','sum r','p','\Deltad','d'); 
axis([0 10 -0.2 1.2]); 
title('Response to a "dividend" shock'); 
print -depsc2 ps2_resp2a.eps; 
