% ps1
% programs for problem set 1

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

addpath('../Data');   % I keep all my data in the same place, and then separate directories for each problem set; 
                      %  this lets matlab find data for load commands. Omit if you're keeping
                      %  data in the same directory as your program
addpath('../Programs');  % Again, I keep programs such as regress_jc all in the same place; omit if you keep them 
                         % in the current directory. 

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

format short; % convenient output formats
format compact; 

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

st = 1;         % 21 for post wwII only; 1 for 26 on. 
plotdata = 1;   % make data plots. 
printdata = 1;  % print data for students 
display = 0;    % display results of each individual regression
createdata = 1; % create data set from raw, or load ps1.txt


% ---------------------------
% simulate AR(1)
% ---------------------------

randn('state',3); %reset random number generator so draw stays the same and you can fool with the plot.
                   % This also allows you to cheat and get plots you like
                   % a sin to which I confess (there's a reason it's 3!)
           
Ext = zeros(20,1);
simsize = 10;
xt = zeros(20,simsize); 
sigxt = zeros(20,1); 
Ext(1) = 1;
xt(1,:) = 1*ones(1,simsize); 
sigxt(1) = 0; 
sige = 0.25;
rho = 0.95; 
for t = 2:20; 
    Ext(t) = rho^(t-1)*Ext(1);  % since we start at 1, j = t-1
    sigxt(t) = ((1-rho^(2*(t-1)))/(1-rho^2))^0.5*sige;
    xt(t,:) = rho*xt(t-1,:)+sige*randn(1,simsize); 
end; 

figure;
plot((0:20)',[0;Ext],'-v','LineWidth',2); 
hold on;
plot((0:20)',0*[0;xt(:,1)],'-k');
axis([0 20 -0.1 1.1])
print -depsc2 ps1_AR(1)IR.eps;

figure;
plot((1:20)',[Ext],'-v','LineWidth',2); 
hold on;
plot((1:20)',[Ext+sigxt ],'-r','LineWidth',2); 
plot((1:20)',[Ext-sigxt ],'-r','LineWidth',2); 
plot((1:20)',[xt],'-k');
plot((1:20)',0*[xt(:,1)],'-k');

legend('E(x_t)','E(x_t)+\sigma(x_t)','E(x_t)-\sigma(x_t)','x_t',3); 
print -depsc2 ps1_simulation.eps



%------------------------
% MA(1) problem
% -----------------------


theta = [0.5 0 -0.5]; %both cases at once
e = [ 0 0 0 1 0 0 0 ]';
t = [-1 0 1 2 3 4 ]';
r = e(2:end)*[1 1 1]+e(1:end-1)*theta;
p = cumsum(r);
figure;
subplot(3,1,1);
plot(t,r(:,1),'-ro');
hold on;
plot(t,p(:,1),'-vb');
plot(t,0*p(:,1),'--k');
title('\theta = +0.5'); 
legend('r','p',2);
axis([-inf inf -0.7 1.7]);

subplot(3,1,2);
plot(t,r(:,2),'-ro');
hold on;
plot(t,p(:,2),'-vb');
title('\theta = 0'); 
axis([-inf inf -0.7 1.7]);
plot(t,0*p(:,1),'--k');

subplot(3,1,3);
plot(t,r(:,3),'-ro');
hold on;
plot(t,p(:,3),'-vb');
title('\theta = -0.5'); 
axis([-inf inf -0.7 1.7]);
plot(t,0*p(:,1),'--k');

print -depsc2 ps1_ma1.eps;
    

if createdata == 1; 
    % -------------------------------------------------------------------------
    % Here's where I make your data file from raw data downloaded from WRDS fyi
    % You don't have this raw data, so don't try to run this part of the
    % program by setting createdate=1! If you really want to run this,
    % email me for the data. 
    % 1. load raw data. 
    % -------------------------------------------------------------------------
    
    % 1. Crsp data. both nyse and all; nasdaq is crazy in end of sample. 
      
    load vwr.txt; % load CRSP value weighted return with and without dividends
    cal = vwr(:,1); 
    r = vwr(:,2); 
    rx = vwr(:,3); 
    

    % --------------------------------------------
    % 2. form dividend price ratio
    %  D_t+1      (P_t+1 + D_t+1 )   P_t
    %  ---         ---------------  ------  - 1 
    %  P_t+1            P_t          P_t+1
    % 
    % = R_t+1 / Rx_t+1 - 1
    
    % then form dividend growth 
    % D_t+1 / D_t = (D_t+1/P_t+1)/(D_t / P_t)*(P_t+1/P_t) 
    % -----------------------------------------------------
    

    dp =  (1+r)./(1+rx) - 1; 
    dd = [NaN ; dp(2:end)./dp(1:end-1).*(1+rx(2:end))-1]; 
    
    cal = cal(st:end); 
    dp = dp(st:end); 
    dd = dd(st:end); 
    r = r(st:end); 

    maxdate = 1925+size(r,1); % last date  in sample
    
    % sbbi for treasury return 
    
    sbbi = load('sbbi_annual.txt'); % load CRSP bills and inflation 
    tbret = sbbi(:,2); 
    tbret = tbret(st:end); 
    
    % print data for problem set 1 in nice text files for students. 
    
    % I use this fancy fprintf call; it makes pretty output but takes a long time to write. Use disp() for first tries
    
    if printdata; 
        fid = fopen('../Data/ps1_data.txt','w'); 
        fprintf(fid,'%% CRSP stock data \n');
        fprintf(fid,'%% all returns and dividend growth are net, i.e. a 10 percent return is 0.10. \n'); 
        fprintf(fid,'%% d/p is a fraction, i.e. a 4 percent d/p ratio is 0.040. \n');      
        fprintf(fid,'%% you can load this file into matlab with the statements \n'); 
        fprintf(fid,'%% load ps1.txt; dates = ps1(:,1); ret = ps1(:,2); dp = ps1(:,3), etc. \n'); 
        fprintf(fid,'%10s %13s %13s %13s %13s \n','% Date','Return','D/P','D growth','T-bill ret.');
        fprintf(fid,'%10i %13.8f %13.8f %13.8f %13.8f \n', [cal r dp dd tbret]'); 
        fclose(fid);
    end;

    % -----------------------------------------------------
    % end raw data construction. All you would have to do 
    % is load the new data set as follows 
    % -----------------------------------------------------
    
    % make extra figure of prices and dividends. allrx isn't saved and not
    % in problem set, so here. 
    
    figure; 
    div = cumsum(log(1+dd(2:end))); 
    prc = cumsum(log(1+rx(2:end))); 
    plot((1927:maxdate)',[div prc]);
    legend('log d','log p',4); 
    axis([1927 maxdate+1 -inf inf]); 
    print -depsc2 ps1b.eps; 

    
end;   % ends if createdata


% now this is how you would do it 
ps1_data = load('..\data\ps1_data.txt'); 
cal = ps1_data(:,1); 
r = ps1_data(:,2); 
dp = ps1_data(:,3); 
dd = ps1_data(:,4); 
tbret = ps1_data(:,5); 

maxdate = 1925+size(r,1); % last date  in sample
disp('sample')
disp([cal(1) cal(end)])

% -------------------
% plots of data
% -------------------

figure; 
subplot(2,1,1); 
plot((1926:maxdate)',dp*100); 
title('Dividend price ratio'); 
axis([1926-1 maxdate+1 0 8]); 
subplot(2,1,2); 
plot((1926:maxdate)',[r*100 tbret*100]); 
title('Returns');
legend('Stock','T Bill',4); 
axis([1926-1 maxdate+1 -50 50]); 
print -depsc2 ps1a.eps; 
% I historically used eps format to import graphs into other applications.
% tiff usually works better, so use print -dtiff ps1a.tif

disp('mean and sd of DP and log DP')
disp([mean(dp) mean(log(dp))]);
disp([std(dp) std(log(dp))]);

disp('mean and sd of returns and dividend growth (logs)'); 
disp([mean(r) mean(dd(2:end))]);
disp([std(r) std(dd(2:end))]); 

% just returns and treasury for lecture
figure; 
subplot(2,1,1); 
plot((1926:maxdate)',[tbret*100 r*100],'linewidth',2); 
title('Stock and T bill return'); 
axis([1926-1 maxdate+1 -inf inf]); 
subplot(2,1,2); 
plot((1926:maxdate)',[tbret*100], 'linewidth',2); 
title('T bill returns expanded');
axis([1926-1 maxdate+1 -inf inf]); 
print -depsc2 ps1_stock_and_tb.eps; 


 figure; 
 plot((1926:maxdate)',[1./dp]); 
 axis([1926-1 maxdate+1 0 inf]);
 title('Price-dividend ratios'); 
 print -depsc2 ps1xa.eps; 
  
figure; 
plot((1926:maxdate)',dp*100,'LineWidth',2); 
title('Dividend price ratio'); 
axis([1926-1 maxdate+1 0 8]); 
print -depsc2 DP_graph.eps; 



% -------------------------------------------
% regressions 
% -------------------------------------------

% 1. regressions of returns on lagged returns

disp('Regression of returns on lagged returns 1926-today'); 
fprintf('       %10s %10s %10s %10s %10s \n','b','t(b)','R2','E(R)','s(E(R))'); 

T = rows(r); 

rhv = [ones(T-1,1) r(1:T-1)]; 
lhv = r(2:end); 
%disp('NYSE return on lagged return, 1926-today'); 
[b,bse,bt,r2,F] = regress_jc(lhv,rhv,display); 
fprintf('Stock  %10.2f %10.2f %10.2f %10.2f %10.2f \n', [b(2) bt(2) r2 100*mean(lhv) 100*std(rhv*b)]);

rhv = [ones(T-1,1) tbret(1:T-1)]; 
lhv = tbret(2:end); 
%disp('T bill return on T bill return, 1926-today'); 
[b,bse,bt,r2,F] = regress_jc(lhv,rhv,display); 
fprintf('T Bill %10.2f %10.2f %10.2f %10.2f %10.2f \n', [b(2) bt(2) r2 100*mean(lhv) 100*std(rhv*b)]);

re = r-tbret; 
rhv = [ones(T-1,1) re(1:T-1)]; 
lhv = re(2:end); 
%disp('NYSE excess return on its lag, 1926-today'); 
[b,bse,bt,r2,F]  = regress_jc(lhv,rhv,display); 
fprintf('Excess %10.2f %10.2f %10.2f %10.2f %10.2f \n', [b(2) bt(2) r2 100*mean(lhv) 100*std(rhv*b)]);

% 2. regressions of excess returns on dividend price ratios

disp(''); 
disp('Regression of excess returns on lagged D/P 1926-today'); 
fprintf('                  %10s %10s %10s %10s %10s %10s %10s \n','b','t(b)','R2','E(R)','s(E(R))','E(R)/k','s(E(R))/k'); 

rhv = [ones(T-1,1) dp(1:T-1)]; 
lhv = re(2:end); 
%disp('NYSE excess return on NYSE D/P ratio, 1926-today'); 
[b,bse,bt,r2,F]  = regress_jc(lhv,rhv,display); 
fprintf('1 Yr.             %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f\n',...
    [b(2) bt(2) r2 100*mean(lhv) 100*std(rhv*b) 100*mean(lhv)/1 100*std(rhv*b)/1]);

r5 = (1+r(1:T-4)).*(1+r(2:T-3)).*(1+r(3:T-2)).*(1+r(4:T-1)).*(1+r(5:T)); 
tbret5 = (1+tbret(1:T-4)).*(1+tbret(2:T-3)).*(1+tbret(3:T-2)).*(1+tbret(4:T-1)).*(1+tbret(5:T)); 
re5 = r5 - tbret5; 
rhv = [ones(T-5,1) dp(1:T-5)]; 
lhv = re5(2:end); % be careful with indices here!
%disp('NYSE 5 year excess return on NYSE D/P ratio, overlapping, 1926-today'); 
[b,bse,bt,r2,F]  = regress_jc(lhv,rhv,display); 
fprintf('5 Yr. overlap     %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f\n',...
    [b(2) bt(2) r2 100*mean(lhv) 100*std(rhv*b) 100*mean(lhv)/5 100*std(rhv*b)/5]);

re5a = re5(2:5:rows(re5)); 
dp5a = dp(1:5:(rows(dp)-5)); 
rhv = [ones(rows(dp5a),1) dp5a]; 
lhv = re5a; 
%disp('NYSE 5 year excess return on NYSE D/P ratio, not overlapping, 1926-today'); 
[b,bse,bt,r2,F]  = regress_jc(lhv,rhv,display); 
fprintf('5 Yr. no overlap  %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f\n',...
    [b(2) bt(2) r2 100*mean(lhv) 100*std(rhv*b) 100*mean(lhv)/5 100*std(rhv*b)/5]);

% doing it right: overlapping data with corrected standard errors -- requires 
% that you get the olsgmm program from the class website. 

lhv1 = (1+r(2:end));
lhv2 = (1+tbret(2:end)); 
rhv = [ones(T-1,1) dp(1:end-1)];
for i = 1:10;
     lhv = lhv1-lhv2; 
     lags = i; 
     weight = 0; 
     [bv(:,i),sebv(:,i),R2v(:,i),R2vadj,v,F] = olsgmm(lhv,rhv,lags,weight);
     % if you got an error message here, go to the class website and get olsgmm. Or just get rid
     % of this section -- this is beyond what I asked you to do. 
     if i == 5; 
         fcst5 = rhv*bv(:,i); 
         actual5 = lhv; 
     end; 
     if i == 10; 
         fcst10 = rhv*bv(:,i); 
         actual10 = lhv;
     end; 
         
     lhv1 = lhv1(1:end-1).*(1+r(2+i:end)); 
     lhv2 = lhv2(1:end-1).*(1+tbret(2+i:end)); 
     rhv = rhv(1:end-1,:); 
 end; 
 
 disp(''); 
 disp('NYSE excess returns on D/P, overlapping, many horizons, 1926-today'); 
 disp('Horizon (years)'); 
 fprintf('   ');  fprintf('%5i ',1:10); fprintf('\n'); 
 fprintf('b  ');  fprintf('%5.2f ',bv(2,:));fprintf('\n'); 
 fprintf('t  ');  fprintf('%5.2f ',bv(2,:)./sebv(2,:));fprintf('\n'); 
 fprintf('R2 ');  fprintf('%5.2f ',R2v);fprintf('\n'); 

 figure; 
 plot((1:10)',bv(2,:),'-rv'); 
 hold on; 
 plot((1:10)',100*R2v,'-bo'); 
 hold on; 
 plot((1:10)',10*bv(2,:)./sebv(2,:),'--g+');
 legend('coeff','R^2 (%)','10 \times t stat',2); 
 xlabel('horizon, years'); 
 print -depsc2 ps1_longhorizonb.eps;
 
 
 plotdates = (1926:maxdate)'; 
 figure; 
 subplot(2,1,1); 
 plot(plotdates(1+5:T),fcst5,'-v', plotdates(1+5:T),actual5,'-');
 title('Actual and forecast 5 year excess returns'); 
% legend('forecast','actual'); 
 axis([1930 maxdate -inf inf]); 
 subplot(2,1,2); 
 plot(plotdates(1+10:T),fcst10,'-v',plotdates(1+10:T), actual10,'-');
 axis([1930 maxdate+1 -inf inf]);
 title('Actual and forecast 10 year excess returns'); 
% legend('forecast','actual');  
 print -depsc2 ps1c.eps; 
 
 
