% 
% Written by:
% -- 
% John L. Weatherwax                2009-04-21
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

close all; drawnow; clear functions; 
rehash; 
clc; 

X = [ 1, -5, 0 ; 
      1, -3, 0 ;  
      1, -1, 0 ; 
      1,  1, 1 ; 
      1,  3, 1 ; 
      1,  5, 1 ; 
];

y = [ 5; 7; 10; 14; 17; 18 ]; 

n = size(X,1); p = 1; % <- note p does not change with the addition of these indication functions

XTX = X' * X; 

% solve the least-squares problem to estimate betas:
fprintf( 'estimated betas are\n' );
beta_hat = XTX \ ( X' * y ) 

% compute the predictions: 
y_hat = X * beta_hat 

% compute the SSE: 
SSE = sum( ( y - y_hat ).^2 ); 

% compute an estimate of \sigma: 
s2 = SSE/(n-p-1); s = sqrt(s2); 

fprintf('covariance matrix of coeficients beta\n');
V_hat_beta_hat = s2 * inv( XTX ) 

% compute an estimate of the standard error of each beta: 
s_beta = sqrt( diag( V_hat_beta_hat ) );

% compute the t statistics (against the hypotheis H_0: \beta_i = 0):
fprintf('t statistics for these coefficients\n');
t_stats = beta_hat ./ s_beta

fprintf('t value for significance\n');
alpha = 0.05; 
tinv( 1 - 0.5*alpha, n-p-1 )


%--
% Refit our regression model without a level change ... we conclude that there is no level change
%--
X = [ 1, -5, ;
      1, -3, ; 
      1, -1, ;
      1,  1, ;
      1,  3, ;
      1,  5, ;
];

y = [ 5; 7; 10; 14; 17; 18 ]; 

n = size(X,1); p = 1; % <- note p does not change with the addition of these indication functions

XTX = X' * X; 

% solve the least-squares problem to estimate betas:
fprintf( 'estimated betas are\n' );
beta_hat = XTX \ ( X' * y ) 

% compute the predictions: 
y_hat = X * beta_hat 

% compute the SSE: 
SSE = sum( ( y - y_hat ).^2 ); 

% compute an estimate of \sigma: 
s2 = SSE/(n-p-1); s = sqrt(s2); 

fprintf('covariance matrix of coeficients beta\n');
V_hat_beta_hat = s2 * inv( XTX ) 

% compute an estimate of the standard error of each beta: 
s_beta = sqrt( diag( V_hat_beta_hat ) );

% compute the t statistics (against the hypotheis H_0: \beta_i = 0):
fprintf('t statistics for these coefficients\n');
t_stats = beta_hat ./ s_beta

fprintf('t value for significance\n');
alpha = 0.05; 
tinv( 1 - 0.5*alpha, n-p-1 )


% lets compute the ANOVA table using this regression
y_bar = mean(y); 
SSTO  = sum( ( y - y_bar ).^2 ); 
SSR   = sum( ( y_hat - y_bar ).^2 ); 
SSE   = SSTO - SSR; 
R2    = SSR/SSTO;

MSR   = SSR / p; 
MSE   = SSE / (n-p-1); 
F     = MSR / MSE; 

finv( 1-alpha, p, n-p-1 )