%%%%%%%%%%%%(c) Jan Kabatek, Universiteit van Tilburg %%%%%%%%%%%%%%%%%%%%%
% This is an auxilliary m-file for 4th tutorial of Econometrics 1 for     %
% research masters.                                                       %
%                                                                         %
% This code is meant to show you how pronounced can be differences        %
% between computation times required by commands used for the same        % 
% mathematical operation                                                               %
%                                                                         %
% IMPORTANT: before asking me about a specific command in this code(or any%
% code) consult a HELP FILE by parsing:	 help <that command into the      %     
% Matlab console.                                                         %	
%																		  %	
% Version 14.10.2013, j.kabatek@uvt.nl, K606 (should you find any typos,  %
% let me know)                                                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
clear all
%%%% data generating process %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 100 ;                           %number of obs
K = 4 ;                             %number of regressors
beta = [1 2 3 4]' ;                 %true regressor values
sigma = .5 ;                        %sd of error terms
 
X = [ones(N,1) randn(N,K-1)] ;      %first col is constant term
E = sigma*randn(N,1) ; 
Y = X*beta+E ; 
%1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nloop = 100000;

display '*** OLS Betas ***'
display 'computation time differences with different inverse functions'
tic
display ' inv(X'''*X)*(X'''*Y)'
for i=1:nloop
    betahat = inv(X'*X)*(X'*Y) ; 
end
toc
display '---------------------------'

tic
display ' (X'''*X)\(X'''*Y) '
for i=1:nloop
    betahat2 = (X'*X)\(X'*Y) ;
end
toc
display '---------------------------'

tic
display ' (Y'''*X)/(X'''*X) ' 
for i=1:nloop
    betahat3 = (Y'*X)/(X'*X)  ;
end
toc
display '---------------------------'

tic
display ' ((X'''*X)^-1)*(X'''*Y) ' 
for i=1:nloop
    betahat4 = ((X'*X)^-1)*(X'*Y) ;
end
toc
display '---------------------------'

tic
display ' (X'''*X)\X'''*Y' 
for i=1:nloop
    betahat5 = (X'*X)\X'*Y ;
end
toc
display '---------------------------'
display '"\" is somewhat faster than the old inv() command, but more importantly - more precise! '
display 'Notice also that the version with bracketed X*Y is much faster than the non-bracketed version (!). '
display 'The intuition for this is as follows:  '
display 'in the former case, we compute inverse of (k x k) matrix X`X, compute (k x 1) matrix X`Y and multiply.'
display 'in the latter case, we take inverse of (k x k) matrix X`X, multiply by (k x n) matrix X` and then by (n x 1) vector Y '
display 'Clearly, the second method performs many more inner multiplications to get to the same result. '
display '---------------------------'
display ' '
display ' '
display ' '
% fyi:
% press ctrl+c to stop the script from running (if there's infinite loop, etc.)
%2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

E = Y - X*((X'*X)\(X'*Y)); %declare E as residuals instead of error terms

display '*** Variance of residuals ***'
display 'computation time differences with different variance functions'
display ' var(E) '
tic 
for i=1:nloop
    var(E);
end
toc
display '---------------------------'
display ' sum(E.^2) '
tic 
for i=1:nloop
    sum(E.^2);
end
toc
display '---------------------------'
display ' E''*E '
tic 
for i=1:nloop
    E'*E;
end
toc
display '---------------------------'
display 'inner product is 5 times faster than the usual var() command! '
display '-> moral of the story is, vectorize and use matrix functions. '
display '---------------------------'

