%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is practice script 1 for the Matlab tutorial%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rand('state',37); % set arbitrary seed for uniform draws
randn('state',37); % set arbitrary seed for normal draws

tic;  % start stop watch

[fid]=fopen('c:\mlab_tut\mlab\logs\script1.txt','w');
if fid==-1;
   warning('File could not be opened');
   break
else;
   disp('File opened successfully');
end;


% generate data
%%%%%%%%%%%%%%%%%%%%%%%%%
n=1000; % set sample size
x1=ones(n,1);
x2= -1.4+randn(n,1);
x3= 3+2*randn(n,1);

X=[x1 x2 x3];
k=size(X,2);
b=[1.2 0.4 0.8]'; 

eps=1.2*randn(n,1);
y=X*b+eps;

bols=inv(X'*X)*X'*y;% get coefficient estimates; inv denotes "inverse"
e=y-X*bols;% Get residuals. Note the asterisk to perform matrix multiplication. 
% Matrix dimensions must be conformable
s2=e'*e/(n-k); %get the regression error (estimated variance of "eps").
Vb=s2*inv(X'*X); % get the estimated variance-covariance matrix of bols
se=sqrt(diag(Vb));% get the standard erros for your coefficients; 
% note the nested command structure: sqrt( ) takes the suare root of
% whatever is in ( ).  Diag ( ) extracts the diagonal from a square matrix
% in ( ).  We can easily combine the two commands.
t=bols./se; % get your t-values.  The dot-operator performs multiplication 
% or division element-by-element.

out=[bols se t]; %combine bols, se, and t vectors into a single matrix
fprintf(fid,'Output table for betas \n'); %label output; the "\n"
% combo is equivalent to "Enter" on your keypad, i.e it moves the next
% piece of output to a new line;
fprintf(fid,'coeff\t\tstd\t\tt-value\n'); % label each column of your output;
% \t inserts extra tabs. Play & experiment with /n and /t until your output 
% looks the way you like it.
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\n',out'); %define the numerical format for
% each column. Here: up to six digits total, with 3 decimal positions. The
% matrix to be plotted (here: "out")appears at the end in TRANSPOSED form
fprintf(fid,'\n'); % this enteres a blank line in your log file

fprintf(fid,'Squared regression error=\t%6.3f \n',s2); % for scalars you can
% define all output formats in a single line
fprintf(fid,'\n');

save c:\mlab_tut\mlab\worksp\script1_stuff y X bols;

finish = toc; %this will show run time in seconds
fprintf(fid,'Time elapsed in seconds \n\n');
fprintf(fid,'%6.3f\n',finish);

st=fclose(fid);
if st==0;
   disp('File closed successfully');
else;
   warning('Problem with closing file');
end;