%---------------------------------------------------
%   OLS and Monte Carlo estimation
%---------------------------------------------------
clear all;
clc;
fid=fopen('output.out','w');
fprintf(fid,'*******************************************************************************\n');
fprintf(fid,'                         output generated by script "ols_mc.m"  \n');
fprintf(fid,'*******************************************************************************\n');
status=fclose(fid);

%----------  # 1 generating data ------------------------
nobs=200;
X=[ones(nobs,1) unifrnd(1,10,nobs,1)];
b=[2;4];
e=2*randn(nobs,1);
y=X*b+e;

%----------  # 2 OLS estimation -------------------------
k=size(X,2);                                % get number of regressors == columns in X
b=inv(X'*X)*X'*y;                           % ols estimator
resid=y-X*b;                                % residuals from OLS
sigma_sq=resid'*resid/(nobs-k);             % variance of residuals
SSE=resid'*resid;                           % Sum of Squared Errors   
SST=(y-mean(y))'*(y-mean(y));               % Sum of Squares Total
R2=1-SSE/SST;                               % R-squared
R2_adj = 1 - (1-R2)*(nobs-1)/(nobs-k);      % R-squared adjusted
var_b=sigma_sq*inv(X'*X);                   % var-covar matrix of coefficients
se_b=sqrt(diag(var_b));                     % s.e.(b)
tstat=b./se_b;                              % t-statistics  
pval=2*(1-tcdf(abs(tstat),nobs-k));         % p-value
CI.a=b-se_b;                                % confidence interval
CI.b=b+se_b;

%---------- printing OLS regression output -------------------------
fid=fopen('output.out','a');
names={'const ' 'X2    '};
fprintf(fid,'\n');
fprintf(fid,'*******************************************************************************\n');
fprintf(fid,'                         Part A: OLS estimation \n');
fprintf(fid,'*******************************************************************************\n');
fprintf(fid,'The equation to be estimated is y = b1*const + b2*X1 + e \n\n');
fprintf(fid,'Number of observations = %6.0f\n',nobs);
fprintf(fid,'Number of parameters to estimate = %5.0f\n',k);
fprintf(fid,'R-squared = %5.5f\n',R2);
fprintf(fid,'R-squared adjusted = %5.5f\n',R2_adj);
fprintf(fid,'SSE = %5.5f\n',SSE);
fprintf(fid,'-------------------------------------------------------------------------------------\n');
fprintf(fid,'Variable    |  Coef.  | Std.Err.   |    t      |  P>|t|   | [95 perc. Conf.Interval]\n');
fprintf(fid,'------------|---------|------------|-----------|----------|--------------------------\n');
for i=1:k
fprintf(fid,'%s', char(names(i)));
fprintf(fid,'      | %5.5f |  %5.5f   |  %5.5f  | %5.5f  | [%5.5f   %5.5f] \n',b(i),se_b(i),tstat(i),pval(i),CI.a(i),CI.b(i));
end
fprintf(fid,'------------------------------------------------------------------------------------\n\n\n');
status=fclose(fid);


%----------  # 3 Monte Carlo estimation -------------------------
% This task can be carried out in two ways:
% either running a loop or using matrix manipulation, which is much faster
% and requires less code.
% The code below exploits both ways, and simulation time for each method
% is recorded into the output file.

%---------- Method 1: Monte Carlo via loop -------------------------
tic;                                % start timer 
M=1000;                             % number of MC runs
b_loop=zeros(M,2);
    for i=1:M
        e=2*randn(nobs,1);
        y=X*b+e;
        b_loop(i,:)=(inv(X'*X)*X'*y)';  % get OLS estimator
    end;
loop.mean_b=mean(b_loop);                 % mean of parameter estimates
loop.sd_b=std(b_loop);                    % standard deviation of parameter estimates
loop.min_b=min(b_loop);                   % minimum of parameter estimates
loop.max_b=min(b_loop);                   % maximum of parameter estimates
toc;                                 % stop timer
t1=toc;                              % record simulation time

%---------- Method 2: Monte Carlo via matrix algebra -------------------------
tic;                        % start timer 
E=2*randn(nobs,M);          % create [nobs x M] matrix of disturbances
Xb=X*b(:,ones(M,1));        % create [nobs x M] matrix in which each column is X*b
Y=Xb+E;
B=(inv(X'*X)*X'*Y)';        % get [M x k] vector of OLS estimators
toc;                        % stop timer
t2=toc;                      % record simulation time

matr.mean_b=mean(B);                  % mean of parameter estimates
matr.sd_b=std(B);                     % standard deviation of parameter estimates
matr.min_b=min(B);                    % minimum of parameter estimates
matr.max_b=min(B);                    % maximum of parameter estimates


%---------- printing Monte Carlo output -------------------------
fid=fopen('output.out','a');
names={'b1    ' 'b2    '};
fprintf(fid,'*******************************************************************************\n');
fprintf(fid,'                  Monte Carlo estimation \n');
fprintf(fid,'*******************************************************************************\n');
fprintf(fid,'\n');
fprintf(fid,'                         Method #1 - using loop \n');
fprintf(fid,'\n');
fprintf(fid,'Number of Monte Carlo simulations:  %6.0f\n',M);
fprintf(fid,'Time required for the simulation:  %5.5f seconds\n',t1);
fprintf(fid,'-----------------------------------------------------------\n');
fprintf(fid,'Variable    |  Mean.  | Std.dev.   |    Min    |  Max     |\n');
fprintf(fid,'------------|---------|------------|-----------|----------|\n');
for i=1:k
fprintf(fid,'%s', char(names(i)));
fprintf(fid,'      | %5.5f |  %5.5f   |  %5.5f  | %5.5f  | \n',loop.mean_b(i),loop.sd_b(i),loop.min_b(i),loop.max_b(i));
end
fprintf(fid,'-----------------------------------------------------------\n\n\n');

fprintf(fid,'\n');
fprintf(fid,'                         Method #2 - using matrices \n');
fprintf(fid,'\n');
fprintf(fid,'Number of Monte Carlo simulations:  %6.0f\n',M);
fprintf(fid,'Time required for the simulation:  %5.5f seconds\n',t2);
fprintf(fid,'-----------------------------------------------------------\n');
fprintf(fid,'Variable    |  Mean.  | Std.dev.   |    Min    |  Max     |\n');
fprintf(fid,'------------|---------|------------|-----------|----------|\n');
for i=1:k
fprintf(fid,'%s', char(names(i)));
fprintf(fid,'      | %5.5f |  %5.5f   |  %5.5f  | %5.5f  | \n',matr.mean_b(i),matr.sd_b(i),matr.min_b(i),matr.max_b(i));
end
fprintf(fid,'-----------------------------------------------------------\n\n\n');
fprintf(fid,'The distributions of parameter estimates are shown in the file output_1.pdf \n');
fprintf(fid,'\n');
status=fclose(fid);

%---------- Plotting histogram of estimates -------------------------
figure(1);
subplot(2,1,1);
hist(b_loop(:,1));
title('Distribution of Constant coefficient estimator: b1');
subplot(2,1,2);
hist(b_loop(:,2));
title('Distribution of X1 coefficient estimator: b2');