%TABLE1   computations for Table 1: 
%         `Financial Asset Returns, 1802-1999'. 
%

%         Ellen McGrattan, 6-15-00

% Subperiods for calculations;

T0       = [1926,1945,1926,1945,1972];
T1       = [1999,1999,1945,1972,1999];
T0       = [1930,1940,1950,1960,1970,1980];
T1       = [1939,1949,1959,1969,1979,1989];

nper     = length(T0);


% Load data

load /home/erm/papers/qr6/data/ibbotson/tablea1.dat
load /home/erm/papers/qr6/data/ibbotson/tablea2.dat
load /home/erm/papers/qr6/data/ibbotson/tablea4.dat
load /home/erm/papers/qr6/data/ibbotson/tablea6.dat
load /home/erm/papers/qr6/data/ibbotson/tablea10.dat
load /home/erm/papers/qr6/data/ibbotson/tablea14.dat
load /home/erm/papers/qr6/data/ibbotson/tablea15.dat
load /home/erm/papers/qr6/data/crsp/wrds/totmkt.dat

% Nominal series

mrlg     = tablea1(:,2:13)';      % 'm' ==> monthly
arlg     = tablea1(:,14);         % 'a' ==> annual,  
mrsm     = tablea4(:,2:13)';      %  if monthly returns in matrix r (nyr x 12)
arsm     = tablea4(:,14);         %  then annual return = (prod(1+r)-1)'
mrgbl    = tablea6(:,2:13)';  
argbl    = tablea6(:,14);
mrgbi    = tablea10(:,2:13)';  
argbi    = tablea10(:,14);
mrgbs    = tablea14(:,2:13)';  
argbs    = tablea14(:,14);
mpi      = tablea15(:,2:13)';
api      = tablea15(:,14);
mrvw     = zeros(12,74);
mrvw(:)  = totmkt(2:889,2);
arvw     = prod(1+mrvw)'-1;


% Real returns (monthly and annual)

mrrsm    = mrsm-mpi;     arrsm  = prod(1+mrrsm)'-1;
mrrlg    = mrlg-mpi;     arrlg  = prod(1+mrrlg)'-1;
mrrvw    = mrvw-mpi;     arrvw  = prod(1+mrrvw)'-1;
mrrgbl   = mrgbl-mpi;    arrgbl = prod(1+mrrgbl)'-1;
mrrgbi   = mrgbi-mpi;    arrgbi = prod(1+mrrgbi)'-1;
mrrgbs   = mrgbs-mpi;    arrgbs = prod(1+mrrgbs)'-1;

% Cumulated Returns 

crsm     = cumprod(1+mrsm(:));
crlg     = cumprod(1+mrlg(:));
crvw     = cumprod(1+mrvw(:));
crgbl    = cumprod(1+mrgbl(:));
crgbi    = cumprod(1+mrgbi(:));
crgbs    = cumprod(1+mrgbs(:));
cpi      = cumprod(1+mpi(:));
crrsm    = cumprod(1+mrrsm(:));
crrlg    = cumprod(1+mrrlg(:));
crrvw    = cumprod(1+mrrvw(:));
crrgbl   = cumprod(1+mrrgbl(:));
crrgbi   = cumprod(1+mrrgbi(:));
crrgbs   = cumprod(1+mrrgbs(:));

nper     = length(T0);
T0       = T0-1925;
T1       = T1-1925;
arith    = zeros(nper,7);
risk     = zeros(nper,7);
value    = zeros(nper,7);
comp     = zeros(nper,7);
rarith   = zeros(nper,6);
rrisk    = zeros(nper,6);
rvalue   = zeros(nper,6);
rcomp    = zeros(nper,6);
for i = 1:nper;

  t0         = T0(i);
  t1         = T1(i);
  lt         = length(T0(i):T1(i));
  A          = [arsm(t0:t1),arlg(t0:t1),arvw(t0:t1), ...
                argbl(t0:t1),argbi(t0:t1),argbs(t0:t1),api(t0:t1)];
  C          = [ones(1,7);crsm,crlg,crvw,crgbl,crgbi,crgbs,cpi];

  Ar         = [arrsm(t0:t1),arrlg(t0:t1),arrvw(t0:t1), ...
                arrgbl(t0:t1),arrgbi(t0:t1),arrgbs(t0:t1)];
  Cr         = [ones(1,6);crrsm,crrlg,crrvw,crrgbl,crrgbi,crrgbs];

  arith(i,:)  = mean(A)*100;
  risk(i,:)   = std(A)*100;
  value(i,:)  = C(12*t1+1,:)./C(12*(t0-1)+1,:);
  comp(i,:)   = (value(i,:).^(1/lt)-1)*100; 

  rarith(i,:) = mean(Ar)*100;
  rrisk(i,:)  = std(Ar)*100;
  rvalue(i,:) = Cr(12*t1+1,:)./Cr(12*(t0-1)+1,:);
  rcomp(i,:)  = (rvalue(i,:).^(1/lt)-1)*100; 

end;

% Add Siegel's estimates:

T0     = T0+1925;
T1     = T1+1925;
T0     = [1802,1802,1871,T0];
T1     = [1997,1870,1925,T1];
nper   = nper+3;

arith  = [NaN, NaN, 9.8, 5.0, NaN, 4.3, 1.3;
          NaN, NaN, 8.1, 4.9, NaN, 5.2, 0.1;
          NaN, NaN, 8.4, 4.4, NaN, 3.8, 0.6;
          arith];
                   
rarith = [NaN, NaN, 8.5, 3.8, NaN, 3.1;
          NaN, NaN, 8.3, 5.1, NaN, 5.4;
          NaN, NaN, 7.9, 3.9, NaN, 3.3;
          rarith];

risk   = [NaN, NaN,17.5, 6.1, NaN, NaN, NaN;
          NaN, NaN,15.5, 2.8, NaN, NaN, NaN;
          NaN, NaN,15.7, 3.0, NaN, NaN, NaN;
          risk];

rrisk  = [NaN, NaN,18.1, 8.8, NaN, 6.1;
          NaN, NaN,16.9, 8.3, NaN, 7.7;
          NaN, NaN,16.8, 6.4, NaN, 4.8
          rrisk];

comp   = [NaN, NaN, 8.4, 4.8, NaN, 4.3, NaN;
          NaN, NaN, 7.1, 4.9, NaN, 5.2, NaN;
          NaN, NaN, 7.2, 4.3, NaN, 3.8, NaN;
          comp];
                   
rcomp  = [NaN, NaN, 7.0, 3.5, NaN, 2.9;
          NaN, NaN, 7.0, 4.8, NaN, 5.1;
          NaN, NaN, 6.6, 3.7, NaN, 3.2;
          rcomp];

value  = [(comp(1:3,:)/100+1).^[[196;69;55]*ones(1,7)];value];
rvalue = [(rcomp(1:3,:)/100+1).^[[196;69;55]*ones(1,6)];rvalue];
 

disp('ANNUAL STOCK MARKET RETURNS, 1802-1999')
disp(' ')

rlab1 = [' ';' '];
rlab2  = [int2str(T0'),ones(nper,1)*'-', ...
          int2str(T1')];
clab2 = ['3    Stocks    ';
         '3 U.S. Treasury';
         '1              '];
clab3 = ['   Small  ';
         '  S&P 500 ';
         ' Value-wgt';
         '  T-Bonds ';
         '  T-Notes ';
         '  T-Bills ';
         ' Inflation'];

disp('NOMINAL ANNUAL RATES OF RETURN')
disp(' ')
disp('Arithmetic average of annual returns:')
disp(mattab1(arith,1,0,0,clab2,clab3,rlab1,rlab2,0))
disp(' ')
disp('Standard deviation of arithmetic returns:')
disp(mattab1(risk,1,0,0,clab2,clab3,rlab1,rlab2,0))
disp(' ')
disp('RESULT OF $1 INVESTED AT START OF PERIOD')
disp('Value at End of Period:')
disp(mattab1(value,1,0,0,clab2,clab3,rlab1,rlab2,0))
disp(' ')
disp('Compounded annual return:')
comp = real(comp);
disp(mattab1(comp,1,0,0,clab2,clab3,rlab1,rlab2,0))
disp(' ')

disp(' ')
disp(' ')

clab2 = ['3    Stocks    ';
         '3 U.S. Treasury'];
clab3 = ['     Small    ';
         '    S&P 500   ';
         '   Value-wgt  ';
         '    T-Bonds   ';
         '    T-Notes   ';
         '    T-Bills   '];

disp('REAL ANNUAL RATES OF RETURN')
disp(' ')
disp('Arithmetic average of annual returns:')
disp(mattab1(rarith,1,0,0,clab2,clab3,rlab1,rlab2,0))
disp(' ')
disp('Standard Deviation of Monthly Return (multiplied by sqrt(12)):')
disp(mattab1(rrisk,1,0,0,clab2,clab3,rlab1,rlab2,0))
disp(' ')
disp('REAL VALUE OF $1 INVESTED AT START OF PERIOD')
disp('Value at End of Period:')
disp(mattab1(rvalue,1,0,0,clab2,clab3,rlab1,rlab2,0))
disp(' ')
disp('Compounded annual return:')
rcomp = real(rcomp);
disp(mattab1(rcomp,1,0,0,clab2,clab3,rlab1,rlab2,0))

