% problem set 4 for asset pricing,

clear all
close all

x = load('..\data\vwr.txt'); 
%x = x(1:end-1,:); recovers weaker dp regression in sample to 2007 -- new
%variable is not affected. 

st = 22; % first date for regressions -- 22 = 194701
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)


% ----------------------------
% 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. 
% ----------------------------

fprintf('start date (first lhv) %8.0f end date %8.0f\n',vw_dates(st), vw_dates(end));


T = size(x,1);
rhv = [ones(T-st+1,1) dp(st-1:T-1)];
lhv = [r(st:T) dd(st:T) dp(st:T)];
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,0,0);
fit1 = rhv*bv; 

fprintf('forecasts using one lag of dp\n') ;
fprintf('   %8.8s %8.8s %8.8s \n', 'b', 't', 'R2')
fprintf('r  %8.3f %8.3f %8.3f \n', bv(2,1),bv(2,1)/sebv(2,1),R2v(1));
fprintf('dd %8.3f %8.3f %8.3f \n', bv(2,2),bv(2,2)/sebv(2,2),R2v(2));
fprintf('dp %8.3f %8.3f %8.3f \n', bv(2,3),bv(2,3)/sebv(2,3),R2v(3));

malags = 5; 
wtv = [0 0.1 0.2 0.3 0.35 0.4 0.5 0.6 0.7]'; 
R2mat = zeros(3,size(wtv,1));
for j = 1:size(wtv);
    wt = wtv(j);
    sumwt = 1; 
    madp = dp(malags:T);
    for i = 2:malags
        madp = madp+wt^(i-1)*dp(malags-i+1:T-i+1);
        sumwt = sumwt + wt^(i-1);
    end;
    madp = madp/sumwt;
    madp = [NaN*ones(malags-1,1) ; madp];

    rhv = [ones(T-st+1,1) dp(st-1:T-1) dp(st-1:T-1)-madp(st-2:T-2)];
    lhv = [r(st:T) dd(st:T) dp(st:T)];
    [bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,0,0);
    R2mat(:,j)=R2v;
end

figure
plot(wtv',R2mat(1:2,:)','-v'); 
legend('return','dd'); 
xlabel('wt'); 
print -depsc2 ps4_r2.eps;

wt = 0.3;
sumwt = 1; 
madp = dp(malags:T);
for i = 2:malags
    madp = madp+wt^(i-1)*dp(malags-i+1:T-i+1);
    sumwt = sumwt + wt^(i-1);
end;
madp = madp/sumwt;
madp = [NaN*ones(malags-1,1) ; madp];

rhv = [ones(T-st+1,1) dp(st-1:T-1) dp(st-1:T-1)-madp(st-2:T-2)];
lhv = [r(st:T) dd(st:T) dp(st:T)];
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,0,0);
fit2 = rhv*bv;

fprintf('forecasts using one lag of dp and %3.0f year MA with weight %3.2f \n',malags,wt) ;
fprintf('   %8.8s %8.8s %8.8s %8.8s %8.8s \n', 'b',  't', 'bma', 'tma', 'R2')
fprintf('r  %8.3f %8.3f %8.3f %8.3f %8.3f \n', bv(2,1),bv(2,1)/sebv(2,1),bv(3,1),bv(3,1)/sebv(3,1),R2v(1));
fprintf('dd %8.3f %8.3f %8.3f %8.3f %8.3f \n', bv(2,2),bv(2,2)/sebv(2,2),bv(3,2),bv(3,2)/sebv(3,2),R2v(2));
fprintf('dp %8.3f %8.3f %8.3f %8.3f %8.3f \n', bv(2,3),bv(2,3)/sebv(2,3),bv(3,3),bv(3,3)/sebv(3,3),R2v(3));

figure
plot([ dp(st-1:T-1)-mean(dp(st-1:T-1)) dp(st-1:T-1)-madp(st-2:T-2) madp(st-2:T-2)-mean(madp(st-2:T-2))])
legend('dp','dp-ma','ma'); 

figure
plot((1926+st:1926+T)',r(st:T),'--k'); 
hold on
plot((1926+st:1926+T)',fit1(:,1),'-r','linewidth',2);
plot((1926+st:1926+T)',fit2(:,1),'-b','linewidth',2);
legend('return',' dp only','dp and ma'); 
title('return forecast')
axis([1945 2010 -inf inf]); 
print -depsc2 ps4_rfcst.eps

figure
plot((1926+st:1926+T)',dd(st:T),'--k'); 
hold on
plot((1926+st:1926+T)',fit1(:,2),'-r','linewidth',2);
plot((1926+st:1926+T)',fit2(:,2),'-b','linewidth',2);
legend('dd',' dp only','dp and ma'); 
title('dd forecast')
axis([1945 2010 -inf inf]); 
print -depsc2 ps4_dfcst.eps

% ----------
% really? comparision with just one lag
% ------------------------

rhv = [ones(T-st+1,1) dp(st-1:T-1) dp(st-1:T-1)-dp(st-2:T-2)];
lhv = [r(st:T) dd(st:T) dp(st:T)];
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,0,0);
fit2 = rhv*bv;

fprintf('forecasts using one lag of dp and dp difference \n') ;
fprintf('   %8.8s %8.8s %8.8s %8.8s %8.8s \n', 'b',  't', 'bD', 'tD', 'R2')
fprintf('r  %8.3f %8.3f %8.3f %8.3f %8.3f \n', bv(2,1),bv(2,1)/sebv(2,1),bv(3,1),bv(3,1)/sebv(3,1),R2v(1));
fprintf('dd %8.3f %8.3f %8.3f %8.3f %8.3f \n', bv(2,2),bv(2,2)/sebv(2,2),bv(3,2),bv(3,2)/sebv(3,2),R2v(2));
fprintf('dp %8.3f %8.3f %8.3f %8.3f %8.3f \n', bv(2,3),bv(2,3)/sebv(2,3),bv(3,3),bv(3,3)/sebv(3,3),R2v(3));

rhv = [ones(T-st+1,1) dp(st-1:T-1) dp(st-1:T-1)-dp(st-2:T-2) dp(st-2:T-2)-madp(st-2:T-2)];
lhv = [r(st:T) dd(st:T) dp(st:T)];
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,0,0);
fit2 = rhv*bv;

fprintf('forecasts using one lag of dp, dp difference, and %3.0f year MA with weight %3.2f \n',malags,wt) ;
fprintf('   %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s \n', 'b',  't', 'b,D',  't,D', 'bma', 'tma', 'R2')
fprintf('r  %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f \n', bv(2,1),bv(2,1)/sebv(2,1),bv(3,1),bv(3,1)/sebv(3,1),bv(4,1),bv(4,1)/sebv(4,1),R2v(1));
fprintf('dd %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f \n', bv(2,2),bv(2,2)/sebv(2,2),bv(3,2),bv(3,2)/sebv(3,2),bv(4,2),bv(4,2)/sebv(4,2),R2v(2));
fprintf('dp %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f \n', bv(2,3),bv(2,3)/sebv(2,3),bv(3,3),bv(3,3)/sebv(3,3),bv(4,2),bv(4,2)/sebv(4,2),R2v(3));



 % ---------------------
 % Goyal Welch Forecast comparision
 % ---------------------

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

x = load('..\data\sbbi.txt'); 
tb_dates = x(:,1); 
Rtb = 1+x(:,3); 

Re = R-Rtb;
T = size(R,1); 

lndp = log(dp);
lnr = log(R); 
 
fcst1 = zeros(T,1);
fcst2 = zeros(T,1);
fcsterr1 = zeros(T,1);
fcsterr2 = zeros(T,1);
coeffs= zeros(T,2); 
for k = 20:T-1;
    b = [ones(k-1,1) lndp(1:k-1)]\lnr(2:k);
    coeffs(k,:)=b';
    fcst1(k+1) = [1 lndp(k)]*b;
    fcsterr1(k+1) = lnr(k+1)-[1 lndp(k)]*b;
    meanr = mean(lnr(1:k));
    fcst2(k+1) = meanr;
    fcsterr2(k+1) = lnr(k+1)-meanr; 
end; 

rmse1 = (cumsum(fcsterr1.^2)./[ones(19,1); (1:T-19)']).^0.5; 
rmse2 = (cumsum(fcsterr2.^2)./[ones(19,1); (1:T-19)']).^0.5; 

disp('root mean squared errors') 
disp('regression, mean'); 
disp(100*[rmse1(T) rmse2(T)]); 

figure;
plot((1927:1:1926+T)',[100*rmse1 100*rmse2]); 
legend('regression','mean',2);
title('root mean squared errors')
print -depsc2 GW_rmse.eps; 

figure;
plot((1927:1:1926+T)',[100*fcsterr1 100*fcsterr2]);     
legend('regression','mean',2);
title('forecast errors'); 

figure; 
plot((1927:1:1926+T)',[100*fcst1 100*fcst2]);     
legend('regression','mean',2);
title('forecasts')
print -depsc2 GW_fcst.eps; 

figure; 
plot((1927:1:1926+T-1)',coeffs(1:T-1,:));     
legend('a','b',2);
title('regression coefficients')    
print -depsc2 GW_coeffs.eps; 


stop

% ------------------------
% A question asking for CAY estimation
% not used 11. Needs update to conform to material moved to lecture notes
% and identities in ir functions.
% maybe compete with lags in dp?  
%-----------------------

close all
clear all

x = load('..\data\all_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

x = load('..\data\sbbi.txt'); 
tb_dates = x(:,1); 
Rtb = 1+x(:,3); 

Re = R-Rtb;
T = size(R,1); 

x = load('..\data\cay_q_07Q4.txt'); 
z = x; 
x = x(4:4:size(x,1),:); % make annual data by choosing end of year value
cay_dates = x(:,1); 
c = x(:,2); 
a = x(:,3); 
y = x(:,4); 
cay = x(:,5); 

% chop vw to same size as cay

vw_dates = vw_dates(27:end);
disp(' '); 
fprintf(' vw data from %8.2f to %8.2f\n',vw_dates(1),vw_dates(end)); 
fprintf('cay data from %8.2f to %8.2f\n',cay_dates(1),cay_dates(end)); 
R = R(27:end); 
dp = dp(27:end); 
Rtb = Rtb(27:end); 
Re = Re(27:end);
DD = DD(27:end); 
T = size(R,1); 

% -------------------------------------------------------
% price and return decompositions
% -------------------------------------------------------

% price and return decomposition using only dp as regressor
rhv = [ones(T-1,1) log(dp(1:T-1)) ]; 
lhv = [log(R(2:T)) log(DD(2:T)) log(dp(2:T)) ];
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,0,0);
err = lhv-rhv*bv;
err_r = [ NaN ; err(:,1)]; 
err_d = [ NaN ; err(:,2)];
err_dp = [ NaN ; err(:,3)];

phi = bv(2,3)'; % VAR coefficients of r dd dp on dp 
br = bv(2,1); 
bd = bv(2,2); 
disp(' '); 
disp('VAR -- dp only '); 
fprintf('    %8.8s %8.8s  \n ','dp', 'R2'); 
fprintf('  r %8.2f %8.2f  \n',bv(2,1),R2v(1)); 
fprintf(' dd %8.2f %8.2f  \n',bv(2,2),R2v(2));
fprintf(' dp %8.2f %8.2f \n',bv(2,3),R2v(3));

% price decomposition 

rho = mean(1./dp)/(1+mean(1./dp));
fprintf('rho %8.3f \n',rho); 

x = log(dp); 
x = x-ones(T,1)*mean(x);
dterm = x*(bd'*inv(eye(1)-rho*phi))'; 
rterm = -x*(br'*inv(eye(1)-rho*phi))';

figure; 
plot((51:1:50+T)',[dterm rterm dterm+rterm -x(:,1)]); 
legend('d','r','d+r','pd',2);
title('var (p-d), dp only') ;
print -depsc2 dp_pdecomp.eps

disp('variance of dp, dp only'); 
components = [dterm rterm -x(:,1)]; 
fprintf('     %8.8s %8.8s %8.8s \n', 'd term','r term','p-d');
fprintf('var  %8.3f %8.3f %8.3f \n',var(components)); 
fprintf('sd   %8.3f %8.3f %8.3f \n',std(components)); 
fprintf('corr \n'); 
fprintf('     %8.3f %8.3f %8.3f \n',corr(components)); 

% return decomposition

Erterm = x*br; 
dterm = rho/(1-rho*phi)*bd*err_dp; 
rterm = -rho/(1-rho*phi)*br*err_dp;

figure; 
plot((1952:1950+T)', ...
    [err_r(2:T) err_d(2:T)...
    dterm(2:T)...
    rterm(2:T)] );
legend('\epsilon_r','\epsilon_d','future d','future r'); 
axis([1950 2010 -inf inf]); 
title('return decomposition, dp only') 
print -depsc2 dp_rdecomp.eps

fprintf('standard deviation and correlation of return components\n'); 
fprintf(' %8.8s %8.8s %8.8s %8.8s %8.8s \n','r','Er','d shock','fut d','fut r'); 
fprintf(' %8.2f %8.2f %8.2f %8.2f %8.2f \n',100*std(log(R(2:end))),100*std(Erterm),100*std(err_d(2:end)),...
                                            100*std(dterm(2:end)),100*std(rterm(2:end))); 
                                        
components = [log(R(2:end)) Erterm(1:end-1) err_d(2:end) dterm(2:end) rterm(2:end)];
fprintf(' %8.2f %8.2f %8.2f %8.2f %8.2f \n',corr(components));
             

% price, dividend and bn trend. 

bp = 1-phi+bd; 
fprintf('bp %8.3f\n',bp); 
dplt = x*bp/(1-phi); 

newdd = log(DD)-mean(log(DD))*ones(T,1); 
newdp = log(dp) - mean(log(dp))*ones(T,1);
newd = [0 ; cumsum(log(DD(2:end)))]; 
newp =  newd - newdp; 

figure; 
plot((1951:1950+T)', [newd newp newp + dplt],'LineWidth',2) 
legend('d','p','BNtrend'); 


% ------------------------------------------------
% now with CAY in the VAR log R and log DD for VAR. 
% -------------------------------------------------

rhv = [ones(T-1,1) log(dp(1:T-1)) cay(1:T-1)]; 
lhv = [log(R(2:T)) log(DD(2:T)) log(dp(2:T)) cay(2:T)];
[bv,sebv,R2v,R2vadj,v,F] = olsgmm(lhv,rhv,0,0);
phi = bv(2:3,3:4)'; % VAR coefficients of r dd dp cay on dp cay 
br = bv(2:3,1); 
bd = bv(2:3,2); 
err = lhv-rhv*bv;
err_r = [ NaN ; err(:,1)]; 
err_d = [ NaN ; err(:,2)];
err_x = [ NaN NaN; err(:,3:4)];


disp(' '); 
disp('VAR'); 
fprintf('    %8.8s %8.8s %8.8s \n ','dp', 'cay', 'R2'); 
fprintf('  r %8.2f %8.2f %8.2f \n',bv(2:3,1),R2v(1)); 
fprintf(' dd %8.2f %8.2f %8.2f \n',bv(2:3,2),R2v(2));
fprintf(' dp %8.2f %8.2f %8.2f \n',bv(2:3,3),R2v(3));
fprintf('cay %8.2f %8.2f %8.2f \n',bv(2:3,4),R2v(4)); 
disp('t') 
fprintf('  r %8.2f %8.2f %8.2f \n',bv(2:3,1)./sebv(2:3,1),R2v(1)); 
fprintf(' dd %8.2f %8.2f %8.2f \n',bv(2:3,2)./sebv(2:3,2),R2v(2));
fprintf(' dp %8.2f %8.2f %8.2f \n',bv(2:3,3)./sebv(2:3,3),R2v(3));
fprintf('cay %8.2f %8.2f %8.2f \n',bv(2:3,4)./sebv(2:3,4),R2v(4)); 


% price variance decomposition in this multivariate VAR

rho = mean(1./dp)/(1+mean(1./dp));
fprintf('rho %8.3f \n',rho); 

x =[log(dp) cay]; 
x = x-ones(T,1)*mean(x);
dterm = x*(bd'*inv(eye(2)-rho*phi))'; 
rterm = -x*(br'*inv(eye(2)-rho*phi))';

disp('phi'); 
disp(phi); 
disp('phi^2'); 
disp(phi^2); 
disp('phi^3'); 
disp(phi^3); 

disp('inv(eye(2) - rho * phi) -- dp and cay var'); 
disp(inv(eye(2)-rho*phi)); 

disp('br''*inv(eye(2)-rho*phi)')
disp(br'*inv(eye(2)-rho*phi))
disp('br'' ')
disp(br')

disp('bd''*inv(eye(2)-rho*phi)')
disp(bd'*inv(eye(2)-rho*phi))
disp('bd'' ')
disp(bd')

figure; 
plot((51:1:50+T)',[dterm rterm dterm+rterm -x(:,1)]); 
legend('d','r','d+r','pd',2); 
title('var (p-d), dp and cay ') ;
print -depsc2 dp_cay_pdecomp.eps

disp('variance of dp, dp and cay '); 
components = [dterm rterm -x(:,1)]; 
fprintf('     %8.8s %8.8s %8.8s \n', 'dterm','rterm','p-d');
fprintf('var  %8.3f %8.3f %8.3f \n',var(components)); 
fprintf('sd   %8.3f %8.3f %8.3f \n',std(components)); 
fprintf('corr \n'); 
fprintf('     %8.3f %8.3f %8.3f \n',corr(components)); 


% return decomposition

Erterm = x*br; 
dterm = rho*err_x*inv(eye(2)-rho*phi)'*bd; 
rterm = -rho*err_x*inv(eye(2)-rho*phi)'*br;

figure; 
plot((1952:1950+T)', ...
    [err_r(2:T) err_d(2:T)...
    dterm(2:T)...
    rterm(2:T)] );
legend('\epsilon_r','\epsilon_d','future d','future r'); 
axis([1950 2010 -inf inf]); 
title('return decomposition, dp & cay') 
print -depsc2 dp_cay_rdecomp.eps

fprintf('standard deviation and correlation of return components\n'); 
fprintf(' %8.8s %8.8s %8.8s %8.8s %8.8s \n','r','Er','d shock','fut d','fut r'); 
fprintf(' %8.2f %8.2f %8.2f %8.2f %8.2f \n',100*std(log(R(2:end))),100*std(Erterm),100*std(err_d(2:end)),...
                                            100*std(dterm(2:end)),100*std(rterm(2:end))); 
                                        
components = [log(R(2:end)) Erterm(1:end-1) err_d(2:end) dterm(2:end) rterm(2:end)];
fprintf(' %8.2f %8.2f %8.2f %8.2f %8.2f \n',corr(components));
       

%  ADD regressions to diagnose the identity
% 1. using VAR long run return and dividend growth coefficients
% 2. and then directly a few leads of dividend growth and returns on both cay
% and dp
