function [ param ] = two_pass(factors,excess_returns)
% Estimate two-pass regression parameters

T=size(excess_returns,1);

first_step_params=([ones(1,T);factors']*[ones(1,T);factors']')\([ones(1,T);factors']*excess_returns);
%first_step_params=inv([ones(1,size(factors',2));factors']*[ones(1,size(factors',2));factors']')*[ones(1,size(factors',2));factors']*excess_returns;

mean_excess_returns=mean(excess_returns)';

h_betas=first_step_params(2:size(first_step_params,1),:);
h_lambdas=(h_betas*h_betas')\(h_betas*mean_excess_returns);
%h_lambdas=inv(h_betas*h_betas')*h_betas*mean_excess_returns;
h_alphas=first_step_params(1,:);
param=[h_alphas';h_betas(:);h_lambdas];
% first come the N alphas, then the N*K betas and lastly the K lambdas

end

