%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), ending period for data (tend),starting period 
%                     for computing moments (t0), 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)
%                     NOTE: use 1,2, ... for dates 
load setdatq
gdp     = gdp/4;
cons    = cons/4;
gc      = gc/4;
gov     = gov/4;
gvtcmp  = gvtcmp/4;
ifixed  = ifixed/4;
iinv    = iinv/4;
netx    = netx/4;
nstate  = 2;      
tbeg    = 1;
tend    = length(gdp);
t0      = 2; 
t1      = tend;
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)
nobs = tend-tbeg+1;


%STEP 2. create additional series if needed
%
y     = gdp-gvtcmp-iinv-netx;
theta = 0.2441;
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]';
[tem1,tem2,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. setup the moment matrix which has dimension nobs x # of moments
%
H    = [kp./y, (kg./y*ones(1,2)).*state,  ...
            (gc./y*ones(1,2)).*state,  ...
            (cons./y*ones(1,2)).*state,  ...
            h, n, h.*h, n.*n];
H    = [H(2:tend,:),(ltech(2:tend)-0.95*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 5. 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 =  Wt/nobs;

%STEP 7. print results in format suitable for Fortran
str=input('Print Ht and Wt? (yes=1,no=0) ');
if str==1;
  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; 
end;

