% ps9_simulation
% uses regression from problem set 1 to illustrate monte carlo/bootstrap

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; 

st = 21; 

load ps1_data.txt; 
cal = ps1_data(st:end,1); 
r = log(1+ps1_data(st:end,2)); 
dp = log(ps1_data(st:end,3)); 
dd = log(1+ps1_data(st:end,4)); 
rf = log(1+ps1_data(st:end,5)); 

disp('sample')
disp([cal(1) cal(end)])

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

% 2. regressions of returns on dividend price ratios

disp(''); 
disp('Regression of log  returns on log  D/P 1926-today'); 
fprintf('                  %10s %10s %10s  \n','b','se(b)','R2'); 
T = size(r,1);
rhv = [ones(T-1,1) dp(1:T-1)]; 
lhv = r(2:end); 
%disp('NYSE return on NYSE D/P ratio, 1926-today');
display = 0; 

[b,bse1,bt,r2,F]  = regress_jc(lhv,rhv,display); 
fprintf('1 Yr. return      %10.3f %10.3f %10.3f \n',...
    [b(2) bse1(2) r2 ]);
[b,bse2,bt,r2,F]  = regress_jc(dp(2:end),rhv,display); 
fprintf('1 Yr. dp          %10.3f %10.3f %10.3f \n',...
    [b(2) bse2(2) r2 ]);
bse = [bse1(2) bse2(2)]';

% redo as vector
b = rhv\[r(2:T) dp(2:T)];
err = [r(2:T) dp(2:T)]-rhv*b; 
disp('correlation matrix of errors'); 
disp(corr(err));
sigma = cov(err);
D = chol(sigma);  % Produces D'*D=sigma hence if e is random N(0,1) then D'*e has covariance matrix D'*D=sigma
if 0; % check
    error_trial = randn(10000,2)*D;
    disp('check errors');
    disp(1000*sigma);
    disp(1000*cov(error_trial));
end

% bootstrap 
Ntrial = 10000; 
bsim = zeros(Ntrial,2); % initialize place to store results
for trial = 1:Ntrial;
    rdpsim = zeros(T,2); % place for simulated data
    rdpsim(1,:) = [r(1) dp(1)]; % initialize with first data point
    selectors = floor(rand(T-1,1)*(T-1)+1); % creates random integers between 1 and T-1, indices for the bootstrap
    error_trial = err(selectors,:);  % pick shocks for this trial by bootstrap, keepng r and dp together
    
    if 0; % used to test -- should get back original series
        for t = 2:T;
            rdpsim(t,:) = [1 rdpsim(t-1,2)]*b+err(t-1,:); % simulate system forward using new errors
        end;
    else
        for t = 2:T;
            rdpsim(t,:) = [1 rdpsim(t-1,2)]*b+error_trial(t-1,:); % simulate system forward using new errors
        end;
    end
    % now regress in simulated data just like actual data
    rhvsim = [ones(T-1,1) rdpsim(1:T-1,2)];
    btemp = rhvsim\rdpsim(2:T,:);
    bsim(trial,:) = btemp(2,:); % only keep slope coefficients
    % you can add here to find the distribution of t statistics, R2, long horizon regressions, you name it.
end
    
figure; 
subplot(2,1,1);
hist(bsim(:,1),50); 
hold on;
plot(b(2,1)*[1 1], [0 500],'-r','linewidth',4);
title('Return coefficient')
subplot(2,1,2);
hist(bsim(:,2),50); 
hold on;
plot(b(2,2)*[1 1], [0 500],'-r','linewidth',4);
title('Dp coefficient')
print -depsc2 ps9_simulation_1.eps

figure;
plot(bsim(:,1),bsim(:,2),'.');
hold on
plot(b(2,1),b(2,2),'sr','markersize',10,'markerfacecolor','r')
title('scatterplot of return and dp coefficients')
xlabel('Return coeffficient')
ylabel('DP coefficient')
print -depsc2 ps9_simulation_2.eps

disp('assumed coefficients')
disp(b(2,:));
disp('mean of simulated coefficients')
disp(mean(bsim));
disp('std errors from regression')
disp(bse')
disp('sd of simulated coefficients -- compare to standard errors')
disp(std(bsim));
disp('percent of sample slope coeffs greater than observed');
disp(sum(bsim > ones(Ntrial,1)*b(2,:))/Ntrial*100);
disp('min and max b observed');
disp(min(bsim))
disp(max(bsim))

% ---------------------
% Monte Carlo
% ---------------------

bsim = zeros(Ntrial,2); % initialize place to store results
for trial = 1:Ntrial;
    rdpsim = zeros(T,2); % place for simulated data
    rdpsim(1,:) = [r(1) dp(1)]; % initialize with first data point
    error_trial = randn(T-1,2)*D; 
    
    for t = 2:T;
            rdpsim(t,:) = [1 rdpsim(t-1,2)]*b+error_trial(t-1,:); % simulate system forward using new errors
    end;
    
    % now regress in simulated data just like actual data
    rhvsim = [ones(T-1,1) rdpsim(1:T-1,2)];
    btemp = rhvsim\rdpsim(2:T,:);
    bsim(trial,:) = btemp(2,:); % only keep slope coefficients
    % you can add here to find the distribution of t statistics, R2, long horizon regressions, you name it.
end
    
figure; 
subplot(2,1,1);
hist(bsim(:,1),50); 
hold on;
plot(b(2,1)*[1 1], [0 500],'-r','linewidth',4);
title('Return coefficient -- monte carlo')
subplot(2,1,2);
hist(bsim(:,2),50); 
hold on;
plot(b(2,2)*[1 1], [0 500],'-r','linewidth',4);
title('Dp coefficient -- monte carlo')
print -depsc2 ps9_monte_1.eps

figure;
plot(bsim(:,1),bsim(:,2),'.');
hold on
plot(b(2,1),b(2,2),'sr','markersize',10,'markerfacecolor','r')
title('scatterplot of return and dp coefficients -- MC')
xlabel('Return coeffficient')
ylabel('DP coefficient')
print -depsc2 ps9_monte_2.eps

disp('');
disp('monte carlo');
disp('mean of simulated coefficients -- monte carlo')
disp(mean(bsim));
disp('std errors from regression')
disp(bse')
disp('sd of simulated coefficients MC -- compare to standard errors')
disp(std(bsim));
disp('percent of sample slope coeffs greater than observed');
disp(sum(bsim > ones(Ntrial,1)*b(2,:))/Ntrial*100);
disp('min and max b observed');
disp(min(bsim))
disp(max(bsim))

% -------------
% bootstrap null of zero 
% -------------

bassumed = [ 0 0.94]';
bsim = zeros(Ntrial,2); % initialize place to store results
for trial = 1:Ntrial;
    rdpsim = zeros(T,2); % place for simulated data
    rdpsim(1,:) = [0 0]; % initialize with first data point
    selectors = floor(rand(T-1,1)*(T-1)+1); % creates random integers between 1 and T-1, indices for the bootstrap
    error_trial = err(selectors,:);  % pick shocks for this trial by bootstrap, keepng r and dp together
    
    for t = 2:T;
        rdpsim(t,:) = rdpsim(t-1,2)*bassumed'+error_trial(t-1,:); % simulate system forward using new errors
    end;
    % now regress in simulated data just like actual data
    rhvsim = [ones(T-1,1) rdpsim(1:T-1,2)];
    btemp = rhvsim\rdpsim(2:T,:);
    bsim(trial,:) = btemp(2,:); % only keep slope coefficients
    % you can add here to find the distribution of t statistics, R2, long horizon regressions, you name it.
end
    
figure; 
subplot(2,1,1);
hist(bsim(:,1),50); 
hold on;
plot(b(2,1)*[1 1], [0 500],'-r','linewidth',4);
title('Return coefficient -- Boot null = 0')
subplot(2,1,2);
hist(bsim(:,2),50); 
hold on;
plot(b(2,2)*[1 1], [0 500],'-r','linewidth',4);
title('Dp coefficient -- Boot null =0')
print -depsc2 ps9_null_1.eps

figure;
plot(bsim(:,1),bsim(:,2),'.');
hold on
plot(b(2,1),b(2,2),'sr','markersize',10,'markerfacecolor','r')
title('scatterplot of return and dp coefficients -- boot null=0')
xlabel('Return coeffficient')
ylabel('DP coefficient')
print -depsc2 ps9_null_2.eps

disp('bootstrap with null = 0'); 
disp('assumed coefficients')
disp(bassumed');
disp('mean of simulated coefficients')
disp(mean(bsim));
disp('std errors from regression')
disp(bse')
disp('sd of simulated coefficients -- compare to standard errors')
disp(std(bsim));
disp('percent of sample slope coeffs greater than observed');
disp(sum(bsim > ones(Ntrial,1)*b(2,:))/Ntrial*100);
disp('min and max b observed');
disp(min(bsim))
disp(max(bsim))