% Program to solve Hands-on 1.1
% written by Marcel Bluhm in January 2012
% prepare machine
clc;
clear all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load data %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[num, txt, raw]= xlsread('DataPWT',3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Estimate Equation (1.2) %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% - form dependent variable
Y=num(:,[1;2;4]);
Y(Y(:,1)==1970,:)=[]; % eliminate first period grgdp because initial
% - form independent variable
X=[num(:,[1;2]),ones(size(num,1),1),num(:,3)];% constant and independent variable
X(X(:,1)==2005,:)=[]; % eliminate last period rgdpl because it could only be used if the
% growth rate from 2005-2010 would be available as
% dependent variable
% - estimate Equation (1.2)
theta=(X(:,3:4)'*X(:,3:4))^(-1)*X(:,3:4)'*Y(:,3);
% - Estimate asymptotic variance: Equation (1.4) and obtain t-stats
u=[Y(:,[1;2]),Y(:,3)-X(:,3:4)*theta]; % error terms
sig_two=1/(size(X,1)-2)*sum(u(:,3).^2); % variance of error terms
Avar_theta_NR=sig_two*(X(:,3:4)'*X(:,3:4))^(-1); % Non-robust variance estimator
t_stat_NR=theta./(diag(Avar_theta_NR).^(.5));
result_NR=[theta,t_stat_NR];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Test for serial correlation in the error terms %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% - Add u to data matrix 'num'
for i=1:size(u,1)
pos=find(num(:,1)==u(i,1) & num(:,2)==u(i,2)); % find right position for current u
num(pos,5)=u(i,3);
end
% - form dependent variable
Y_=num(:,[1;2;4]);
Y_(Y_(:,1)==1970|Y_(:,1)==1975,:)=[]; % eliminate first two periods of grgdp
% - form independent variable
X_=[num(:,[1;2]),ones(size(num,1),1),num(:,[3;5])];% constant and independent variables
X_(X_(:,1)==2005|X_(:,1)==1970,:)=[]; % eliminate first and last period rgdpl
% - estimate Equation (1.2)
theta_=(X_(:,3:5)'*X_(:,3:5))^(-1)*X_(:,3:5)'*Y_(:,3);
% - compute heteroscedasticity robust variance for theta_
u_=Y_(:,3)-X_(:,3:5)*theta_;
Psi_=diag(diag(u_*u_'));
Avar_theta_=(X_(:,3:5)'*X_(:,3:5))^(-1)*X_(:,3:5)'*Psi_*X_(:,3:5)*(X_(:,3:5)'*X_(:,3:5))^(-1);
% - compute t-statistic for the coefficient of lagged u
t_stat_=theta_(3)/sqrt(Avar_theta_(3,3));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Test for heteroscedasticity in the error terms %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% - form dependent variable
Yh=[num(:,1:2),num(:,5).^2];
Yh(Yh(:,1)==1970,:)=[]; % eliminate first period of uit
% - form independent variable
Xh=[num(:,[1;2]),ones(size(num,1),1),num(:,3),num(:,3).^2];% constant and independent variables
Xh(Xh(:,1)==2005,:)=[]; % eliminate last period rgdpl
% - Estimate POLS of Equation (8) and obtain residuals
thetah=(Xh(:,3:5)'*Xh(:,3:5))^(-1)*Xh(:,3:5)'*Yh(:,3);
epsilon_hat=Yh(:,3)-Xh(:,3:5)*thetah; % residuals from estimating Equation (8)
% - Regress Yh on a constant only
thetah=(Xh(:,3)'*Xh(:,3))^(-1)*Xh(:,3)'*Yh(:,3);
epsilon_tilde_hat=Yh(:,3)-Xh(:,3)*thetah; % residuals from estimating Equation (8) with a constant only
% Obtain R2c
R2ch=1-(epsilon_hat'*epsilon_hat)/(epsilon_tilde_hat'*epsilon_tilde_hat);
test_stath=size(Yh,1)*R2ch;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute Robust Standard Errors following Equation (5) %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p1=zeros(2,2); % initialize first part of Equation (5)
p2=zeros(2,2); % initialize second part of Equation (5)
for i=1:num(end,2) % counter over the cross-sections
Xi=X(X(:,2)==i,3:4);
ui=u(u(:,2)==i,3);
p1=p1+Xi'*Xi;
p2=p2+Xi'*(ui*ui')*Xi;
end
Avar_theta=p1^(-1)*p2*p1^(-1);
t_stat=theta./(diag(Avar_theta).^(.5));
result=[theta,t_stat];