% Program to solve Hands-on 1.3
% written by Marcel Bluhm in January 2012
% prepare machine
clc;
clear all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load data %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[num, txt, raw]= xlsread('DataPWT',3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Form FE transformed equation %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% - form dependent variable
Y=num(:,[1;2;4]);
Y(Y(:,1)==1970,:)=[]; % eliminate first period grgdp
% - form independent variable
X=[num(:,[1;2]),num(:,3)];% independent variable
X(X(:,1)==2005,:)=[]; % eliminate last period rgdpl
N=num(end,2); % number of cross-sections
T=sum(numel(X(X(:,2)==1)),1); % number of observations per cross-section
for i=1:N
X(X(:,2)==i,3)=X(X(:,2)==i,3)-mean(X(X(:,2)==i,3));
Y(Y(:,2)==i,3)=Y(Y(:,2)==i,3)-mean(Y(Y(:,2)==i,3));
end
% - estimate Equation (1.3)
beta=(X(:,3)'*X(:,3))^(-1)*X(:,3)'*Y(:,3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Obtain residuals and estimate avar_beta %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u=[Y(:,[1;2]),Y(:,3)-X(:,3)*beta];
sig_u=1/(N*(T-1)-1)*(u(:,3)'*u(:,3));
p1=0; % initialize summation part from Equation (16)
for i=1:N
p1=p1+X(X(:,2)==i,3)'*X(X(:,2)==i,3);
end
avar_beta=sig_u*p1^(-1);
% avar_beta=sig_u*(X(:,3)'*X(:,3))^(-1); % equivalent to above without loop
t_stat=beta./(diag(avar_beta).^(.5));
result=[beta,t_stat];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Test for serial correlation %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% make auxialiary regression of u_iT on u_i,T-1
Ya=u(u(:,1)==2005,3); % T=2005
Xa=u(u(:,1)==2000,3); % T-1=2000
delta=(Xa'*Xa)^(-1)*Xa'*Ya; % coefficient estimate
e=Ya-Xa*delta; % error term of auxiliary regression
sig_e=1/(numel(Ya)-1)*(e'*e) % variance of error terms of auxiliary regression
avar_delta=sig_e*(Xa'*Xa)^(-1); % variance of delta
t_stat=(delta+1/(T-1))/sqrt(avar_delta);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute robust variance estimator of beta %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p2=0; % initialize second part of Equation (1.17)
for i=1:N
Xi=X(X(:,2)==i,3);
ui=u(u(:,2)==i,3);
p2=p2+Xi'*(ui*ui')*Xi;
end
avar_beta_robust=(X(:,3)'*X(:,3))^(-1)*p2*(X(:,3)'*X(:,3))^(-1);
t_stat=beta/avar_beta_robust^(.5);
result_fe_robust=[beta,t_stat];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Implement FEGLS Procedure %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u(u(:,1)==2005,:)=[]; % drop the last time period
X(X(:,1)==2000,:)=[]; % drop the last time period
Y(Y(:,1)==2005,:)=[]; % drop the last time period
% - obtain consistent estimate for omega
p1=zeros((T-1),(T-1)); % initialize summation part for omega
for i=1:N
u_i=u(u(:,2)==i,3);
p1=p1+u_i*u_i';
end
omega=N^(-1)*p1;
% - estimate Equation (18)
p1=zeros((size(X,2)-2),(size(X,2)-2)); % initialize part one of Equation (18)
p2=zeros((size(X,2)-2),1); % initialize part two of Equation (18)
for i=1:N
Xi=X(X(:,2)==i,3);
Yi=Y(Y(:,2)==i,3);
p1=p1+Xi'*omega^(-1)*Xi;
p2=p2+Xi'*omega^(-1)*Yi;
end
beta_fegls=p1^(-1)*p2;
avar_beta_fegls=p1^(-1);
t_stat=beta_fegls/avar_beta_fegls^(.5);
result_fegls=[beta_fegls,t_stat];