%% Frisch-Waugh theorem (Michele Piffer, LSE, 23/11/2011)
clc 
clear 

%% Generate DGP

n = 100;                % number of rows
k1 = 5;                 % number of first regressors
X1 = 10*rand(n,k1);     % random values assigned
k2 = 3;                 % number of second regressors
X2 = 10*rand(n,k2);     % random values assigned
X = [X1 X2];            % merge regressors in a unique X matrix

epsilon = randn(n,1);   % true errors

beta_true = [4 7 11 23 -8 -4 23 -11]'; % true parameters, which you want to estimate

y = X*beta_true + epsilon; 

%% Estimate model using X and y

beta_hat = (X'*X)^(-1)*X'*y;    % OLS regressors using full X

display(beta_true')
display(beta_hat')     % compare values, note that they are  ...
                       ...quite close to true values

%% Estimate betas on X1 using FW theorem

% regress X1 on X2, get residuals
step1 = (eye(n) - X2*(X2'*X2)^(-1)*X2')*X1;

% regress y on X2, get residuals
step2 = (eye(n) - X2*(X2'*X2)^(-1)*X2')*y;

% regress y on step1
beta_hat_FW1 = (step1'*step1)^(-1)*step1'*y; 

% regress step2 on step1
beta_hat_FW2 = (step1'*step1)^(-1)*step1'*step2;

display(beta_true')
display(beta_hat') 
display(beta_hat_FW1')
display(beta_hat_FW2') % of course, the same on X2 and no
                  