% ************************************************************************ 
% HALL'S RANDOM-WALK HYPOTHESIS -- COMPARISON OF OLS, 2SLS & GMM ESTIMATES 
% ************************************************************************ 

clc;
global n dc dc1 dc2 dc3 dc4 dy;

% Load Data 
M = importdata('matlabex18data.txt');
data = M.data;

% Define Variables 
y = data(1:228,1);
c = data(1:228,5);
n = size(c,1);

% Create Changes in Variables 
dc = c(2:n) - c(1:n-1);
dy = y(2:n) - y(1:n-1);
n = size(dc,1);
xmat = dy;
k = size(xmat,2);

% ************** 
% OLS ESTIMATION 
% ************** 

b = inv(xmat'*xmat)*(xmat'*dc);
resids = dc - xmat*b;
s2 = (resids'*resids)/(n - k);
varb = s2*inv(xmat'*xmat);
tols = b ./ diag(sqrt(varb));

disp(['OLS Est.   tstat']);
disp([[b,tols]]);

% *************** 
% 2SLS ESTIMATION 
% *************** 

% Lagged Changes 
dc1 = dc(4:n-1);
dc2 = dc(3:n-2);
dc3 = dc(2:n-3);
dc4 = dc(1:n-4); 
dc = dc(5:n);
dy = dy(5:n);
n = size(dy,1);

% Estimation 
wmat = [dc1,dc2,dc3,dc4];
dyhat = wmat*inv(wmat'*wmat)*(wmat'*dy);
zmat = dyhat;
betahat = inv(zmat'*zmat)*(zmat'*dc);
resids = dc - dy*betahat;
sigma2hat = (1/n)*(resids'*resids);
varbetahat = sigma2hat*inv(zmat'*zmat);
t2sls = betahat ./ diag(sqrt(varbetahat));

disp(['2SLS Est.  tstat']);
disp([[betahat,t2sls]]);

% ************** 
% GMM ESTIMATION 
% ************** 

% Optimization
bgmm = fminbnd(@GMM,0,1);
	  
% Variance-Covariance Matrix 
m1 = dc1.*(dc - bgmm*dy);
m2 = dc2.*(dc - bgmm*dy);
m3 = dc3.*(dc - bgmm*dy);
m4 = dc4.*(dc - bgmm*dy);
m = [m1,m2,m3,m4];
W = inv((1/n)*(m'*m));

g1 = mean(-dc1.*dy);
g2 = mean(-dc2.*dy);
g3 = mean(-dc3.*dy);
g4 = mean(-dc4.*dy);
capgamma = [g1,g2,g3,g4];
gmmvar = (1/n)*inv(capgamma*W*capgamma');
tgmm = bgmm ./ diag(sqrt(gmmvar));

disp(['GMM Est.   tstat']);
disp([[bgmm,tgmm]]);

% Test of Overidentifying Restrictions 
Wald = n*GMM(bgmm);
disp(['Wald Statistic (df = 3) = ' num2str(Wald')]);