%SETHT    is a script file that calculates moment restrictions needed to
%         implement the procedure in Ingram and Lee.
%         The file returns:
%             (1)  Ht-the average of the moment conditions (1xnmoments)
%             (2)  Wt-the variance covariance matrix of the disturbance 
%                   (nmoments x nmoments)
%                 IMPORTANT NOTE: Wt needs to be scaled by (1+tend/N)
% 

%INITIAL CONDITIONS:  set # of states (nstate), starting period for data
%                     (tbeg, e.g., 1929), ending period for data (tend),
%                     starting period for computing moments (t0, e.g. 
%                     t0=1941), ending period for computing moments (t1),
%                     number of autocorrelations to be used in Wt (ncor), 
%                     switch that determines if autocorrelations are to be 
%                     damped (idamp=1 or 0), 
nstate  = 4;      
tbeg    = 1929;
tend    = 1985;
t0      = 1941; 
t1      = 1985;
ncor    = 3;
idamp   = 1;


%STEP 1. load in data, (e.g., cons,gov,ifixed,iinv,netx,gdp,pop,gvtcmp, ...
%                             n,h,kp,ig,kg,sold,gdef,state)
load setdat
kp=kp+kg;
nobs = tend-tbeg+1;


%STEP 2. create additional series if needed
%
y     = gdp-gvtcmp-iinv-netx;
tech  = y./((kp+kg).^(theta).*n.^(1-theta).*h);
tem   = ols([ones(nobs,1),[1:nobs]'],log(tech));
ltech = log(tech)-tem(1)*ones(nobs,1)-tem(2)*[1:nobs]';
[tem,tem,epst] = ols(ltech(1:nobs-1),ltech(2:nobs));


%STEP 3. if nstate=5, reset state matrix
%
if nstate==5;
  state=[state,zeros(nobs,1)];
  state(13:14,:)=[zeros(2,4),[1;1]];
end


%STEP 4. set up some parameters (not changed)
tend = tend-tbeg+1;
t0   = t0-tbeg+1;
t1   = t1-tbeg+1;
tbeg = 1;


%STEP 5. setup the moment matrix which has dimension nobs x # of moments
%(kg./y*ones(1,4)).*state(:,2:5),  ...
H    = [kp./y,  ...
            (gc./y*ones(1,3)).*state(:,2:4),  ...
            (h*ones(1,3)).*state(:,2:4),  ...
            cons./y, n, h.*h, n.*n];
H    = [H(2:tend,:),(ltech(2:tend)-0.81*ltech(1:tend-1)).^2, ...
            ltech(2:tend).*ltech(1:tend-1)];
lost = 1;
nobs = t1-t0+1;
H    = H([t0:t1]-lost,:);
Ht   = sum(H)'/nobs;


%STEP 6. construct Wt using Newey-West

u     = H-ones(nobs,1)*Ht';
Wt    = u'*u;
for ii= 1:ncor;
 uu   = u(ii+1:nobs,:)'*u(1:nobs-ii,:);
 if idamp == 1
   Wt = Wt + (1.0-ii/ncor)*(uu+uu');
 else
   Wt = Wt + uu+uu';
 end
end
Wt =  inv(Wt/nobs);

%STEP 7. print results in format suitable for Fortran
nmom=length(Ht);
for i=1:nmom;
  for j=1:nmom;
    disp(sprintf('      wt(%d,%d) = %g',i,j,Wt(i,j)))
  end;
end;
for i=1:nmom;
    disp(sprintf('      hdat(%d)  = %g',i,Ht(i)))
end; 

