@ INDIVIS.G @    @ This file is for Version 3 of GAUSS-386i @

@ Craig Burnside - December 1992 @

@ opening procedure to solve indivisible labor model @

proc indivis(beta,alpha,lngamm,deltak,rho) ;
local nc,ns,nex,nf,gammax,kyratio,iyratio,cyratio,mu,muk,mun,mcc,mcs,
 mce,mss0,mss1,msc0,msc1,mse0,mse1,fc,fx,fe,w,r,q,p0,lamb0,lambr,pr,pim,
 alamb,lambz,lambs,lambda,p,lamb1,lamb2,p11,p12,p21,p22,ps,ps11,ps12,ps21,
 ps22,rxe,rle,qxe,qle,phi0,phi1,psi,i,xx,xe,solx,lx,lex,soll,cxl,ce,solc,
 solf,m,h ;

nc = 2 ;
ns = 1 ;
nex = 1 ;
nf=4 ;

gammax=exp(lngamm) ;
kyratio=(1-alpha)/(gammax/beta-(1-deltak)) ;
iyratio=(gammax-(1-deltak))*kyratio ;
cyratio=1-iyratio ;

mu=(gammax-beta*(1-deltak))/gammax;
muk=-alpha*mu ;
mun=alpha*mu ;

mcc = ((-1)~(0))|
      (0~(1-alpha)) ;
mcs = (0~1)|
      ((1-alpha)~1) ;
mce = 0|
      1 ;
mss0 = (muk~1)|
       ((-gammax*kyratio)~0) ;
mss1 = (0~(-1))|
       ((1-alpha+(1-deltak)*kyratio)~0) ;
msc0 = (0~(-mun))|
       (0~0) ;
msc1 = (0~0)|
       (cyratio~(-alpha)) ;
mse0 = (-mu)|
       0 ;
mse1 = 0|
       (-1) ;

fc = (0~alpha)|
     (0~(alpha-1))|
     ((-cyratio/iyratio)~(alpha/iyratio))|
     (0~mun) ;
fx = (1-alpha)|
     (1-alpha)|
     ((1-alpha)/iyratio)|
     muk ;
fe = 1|
     1|
     (1/iyratio)|
     mun ;

w = -inv(mss0 - msc0*inv(mcc)*mcs)*(mss1 - msc1*inv(mcc)*mcs);
r =  inv(mss0 - msc0*inv(mcc)*mcs)*(mse0 + msc0*inv(mcc)*mce);
q =  inv(mss0 - msc0*inv(mcc)*mcs)*(mse1 + msc1*inv(mcc)*mce);
{lamb0,p0} = eigv(w);
if sumc(abs(imag(lamb0))) gt 1e-10 ;
 print "Error - some eigenvalues of W are complex numbers" ;
 print "Results will be misleading" ;
endif ;
pr=real(p0)./(sqrt(sumc(real(p0)^2))') ;
lambr=real(lamb0) ;
alamb=abs(lambr) ;

lambz=sortind(alamb) ;
lambs=lambr[lambz,1] ;
lambda=diagrv(eye(2*ns),lambs) ;
p=pr[.,lambz] ;

lamb1 = lambda[1:ns,1:ns] ;
lamb2 = lambda[ns+1:2*ns,ns+1:2*ns] ;

p11 = p[1:ns,1:ns] ;
p12 = p[1:ns,ns+1:2*ns] ;
p21 = p[ns+1:2*ns,1:ns] ;
p22 = p[ns+1:2*ns,ns+1:2*ns] ;

ps=inv(p) ;
ps11 = ps[1:ns,1:ns] ;
ps12 = ps[1:ns,ns+1:2*ns] ;
ps21 = ps[ns+1:2*ns,1:ns] ;
ps22 = ps[ns+1:2*ns,ns+1:2*ns] ;

rxe = r[1:ns,1:nex] ;
rle = r[ns+1:2*ns,1:nex] ;
qxe = q[1:ns,1:nex] ;
qle = q[ns+1:2*ns,1:nex] ;

phi0=ps21*rxe+ps22*rle ;
phi1=ps21*qxe+ps22*qle ;

psi=zeros(ns,nex) ;
i=1 ;
do while i le ns ;
 psi[i,.] = -(phi0[i,.]*rho+phi1[i,.])*inv(eye(nex)-rho./lamb2[i,i])/
              lamb2[i,i] ;
 i=i+1 ;
endo ;

xx = p11*lamb1*inv(p11) ;
xe = (p11*lamb1*ps12+p12*lamb2*ps22)*inv(ps22)*psi+rxe*rho+qxe ;
solx = xx~xe ;

lx = -inv(ps22)*ps21 ;
lex = inv(ps22)*psi ;
soll = lx~lex ;

cxl = inv(mcc)*mcs ;
ce = inv(mcc)*mce ;
solc = ( cxl*(eye(ns)|lx) )~( cxl*(zeros(ns,ns)|lex)+ce ) ;

solf = (fx~fe)+fc*solc ;

m =  solx|(zeros(nex,ns)~rho) ;
h = soll|solc|solf ;
retp(m|h) ;
endp ;

proc hpmom(sigma,drules,ncorr) ;
local nex,ns,m,h,sigh,dr,di,vr,vi,sigt,gammt0,gamm0,hbig,base,maxf,jp,jm,
 hpap,hpam,hpa,ny,gammf,facv,i,j,gammj,tcorr,facr,sd,vr1,dr1 ;
nex=rows(sigma) ;
ns=cols(drules)-nex ;
m=drules[1:nex+ns,.] ;
h=drules[nex+ns+1:rows(drules),.] ;
sigh=zeros(ns+nex,ns+nex) ;
sigh[ns+1:ns+nex,ns+1:ns+nex]=sigma ;

{dr1,vr1}=eigv(m) ;
if sumc(abs(imag(dr1))) gt 1e-10 ;
 print "Error - some eigenvalues of M are complex numbers" ;
 print "Results will be misleading" ;
endif ;
vr=real(vr1)./(sqrt(sumc(real(vr1)^2))') ;
dr=real(dr1) ;

sigt=inv(vr)*sigh*inv(vr') ;
gammt0=(ones(ns+nex,ns+nex)./(ones(ns+nex,ns+nex)-dr.*(dr'))).*sigt ;
gamm0=vr*gammt0*(vr') ;
hbig=eye(ns+nex)|h ;

base=181;
maxf=101;

jp=seqa(1,1,maxf) ;
jm=seqa(maxf,-1,maxf) ;
hpap=-(.894^jp).*(.0561*cos(.112*jp)+.0558*sin(.112*jp)) ;
hpam=-(.894^jm).*(.0561*cos(.112*jm)+.0558*sin(.112*jm)) ;
hpa=hpam|(1-(.0561*cos(0)+.0558*sin(0)))|hpap ;

ny=rows(drules) ;
gammf=zeros(1,ncorr+1).*.zeros(ns+nex,ns+nex) ;
facv=zeros(ny,(ncorr+1)*ny) ;
i=0 ;
do while i le ncorr ;
 j=0 ;
 do while j le base ;
  if j eq 0 ;
   gammj=gamm0 ;
   gammf[.,i*(ns+nex)+1:(i+1)*(ns+nex)]=gammf[.,i*(ns+nex)+1:(i+1)*(ns+nex)]+
    gammj*(hpa[i+1:2*maxf+1,1]'*hpa[1:2*maxf+1-i,1]) ;

  else ;
   gammj=m*gammj ;
   if j le i ;
    gammf[.,i*(ns+nex)+1:(i+1)*(ns+nex)]=gammf[.,i*(ns+nex)+1:(i+1)*(ns+nex)]+
        gammj*(hpa[i+1-j:2*maxf+1,1]'*hpa[1:2*maxf+1-i+j,1])+
        (gammj')*(hpa[i+j+1:2*maxf+1,1]'*hpa[1:2*maxf+1-i-j,1]) ;
   else ;
    gammf[.,i*(ns+nex)+1:(i+1)*(ns+nex)]=gammf[.,i*(ns+nex)+1:(i+1)*(ns+nex)]+
        gammj*(hpa[1:2*maxf+1+i-j,1]'*hpa[1-i+j:2*maxf+1,1])+
        (gammj')*(hpa[i+j+1:2*maxf+1,1]'*hpa[1:2*maxf+1-i-j,1]) ;
   endif ;
  endif ;
  j=j+1 ;
 endo ;
 facv[.,i*ny+1:(i+1)*ny]=hbig*gammf[.,i*(ns+nex)+1:(i+1)*(ns+nex)]*(hbig') ;
 i=i+1 ;
endo ;
sd=sqrt(diag(facv[1:ny,1:ny])) ;
tcorr=sd.*(sd') ;
facr=facv./(ones(1,ncorr+1).*.tcorr) ;
retp(facv|facr) ;
endp ;

proc hpfilter(y,lamb) ;
local d,a,i ;
d=rows(y) ;
a=zeros(d,d) ;
i=3 ;
do while i le d-2 ;
 a[i,i]=6*lamb+1;
 a[i,i+1]=-4*lamb;
 a[i,i+2]=lamb;
 a[i,i-2]=lamb;
 a[i,i-1]=-4*lamb;
 i=i+1 ;
endo ;
a[2,2]=1+5*lamb;
a[2,3]=-4*lamb;
a[2,4]=lamb;
a[2,1]=-2*lamb;
a[1,1]=1+lamb;
a[1,2]=-2*lamb;
a[1,3]=lamb ;
a[d-1,d-1]=5*lamb+1;
a[d-1,d]=-2*lamb;
a[d-1,d-2]=-4*lamb;
a[d-1,d-3]=lamb;
a[d,d]=1+lamb;
a[d,d-1]=-2*lamb;
a[d,d-2]=lamb;
retp((eye(d)-inv(a))*y) ;
endp ;

proc moments(b) ;
local beta,alpha,lngamm,deltak,rho,sigma,drules,ny,mms,sd,facr,mm ;
beta=1.03^(-.25) ;
alpha=b[8,1] ;
lngamm=b[6,1] ;
deltak=b[7,1] ;
rho=b[3,1] ;
sigma=b[4,1] ; sigma=sigma^2 ;
drules=indivis(beta,alpha,lngamm,deltak,rho) ;
ny=rows(drules) ;
mms=hpmom(sigma,drules,0) ;
sd=sqrt(diag(mms[1:ny,.])) ;
facr=mms[ny+1:2*ny,.] ;
mm=sd[6,1]|(sd[4,1]/sd[6,1])|(sd[8,1]/sd[6,1])|(sd[5,1]/sd[6,1])|
(sd[5,1]/sd[7,1])|facr[5,7] ;
retp(mm) ;
endp ;

proc datamom(b) ;
local dm ;
dm=b[9:13,1]|((b[9,1]^2-(b[12,1]^2*b[9,1]^2)-b[13,1]^(-2)*b[12,1]^2*
   b[9,1]^2)/(2*b[13,1]^(-1)*b[12,1]^2*b[9,1]^2)) ;
retp(dm) ;
endp ;

varb=0 ;

proc tests(b) ;
local nm,mm,gm,dm,gd,vm,vd,testn,testd,test,pv,i ;
mm=moments(b) ;
nm=rows(mm) ;
gm=grad2(&moments,b,mm,0) ;
dm=datamom(b) ;
gd=grad2(&datamom,b,dm,0) ;
load varb ;
vm=zeros(nm,1) ;
vd=zeros(nm,1) ;
test=zeros(nm,1) ;
pv=zeros(nm,1) ;
i=1 ;
do while i le nm ;
 vm[i,1]=gm[i,.]*varb*(gm[i,.]') ;
 vd[i,1]=gd[i,.]*varb*(gd[i,.]') ;
 testn=mm[i,1]-dm[i,1] ;
 testd=(gm[i,.]-gd[i,.])*varb*((gm[i,.]-gd[i,.])') ;
 test[i,1]=(testn^2)/testd ;
 pv[i,1]=cdfchic(test[i,1],1) ;
 i=i+1 ;
endo ;
format /rd /m1 11,4 ;
print "Model Moments   s.e.      Data Moments     s.e       Test      P-Value";
print mm~sqrt(vm)~dm~sqrt(vd)~test~pv ;
retp(0) ;
endp ;

   /* To run this example file, type
       RUN MINQUAD.SET [F4][F2]
       RUN GMMQ.SET [F4][F2]
       RUN INDIVIS.G [F4][F2]  */

@ ======================== User Definition Area ============================ @
@ ********************************************************************** @
@ ----------- Step 1: Prepare Output File ---------------------------@
@ ********************************************************************** @
output file=indivis.out reset; @ Specify the name of the output file. @

? "************************** GMM Results *********************************";
? "INDIVIS.OUT";
? "Indivisible Labor Model Deterministic Trend ";
datestr(0);
timestr(0);
? "Parameters  'theta,lna,rhoa,siga,ay,lngamm,delta,alpha,5 data mom'";

@ *********************************************************************** @
@------------ Step 2:  Define Global Variables  ------------------------- @
@ *********************************************************************** @

@----------------------------@
   @ Section A. User must specify variables in this section. @
@----------------------------@

tend=113;   @ # of observations=T; scalar @

bgm=1|1|0.9|0.01|1|.01|.02|0.6|.01|1|1|1|1 ;
                @ bgm=b(1)|..|b(kgm); kgm by 1 vector
                      Initial values of the coefficients
                 (In this example, the parameters happen to converge to
                  values that are not economically meaningful when W0=I.
                  The initial values are set to values close to the converged
                  valuse in this example, so that this example converges
                  quickly.) @

nw=13;           @ # of disturbance terms in w(t); scalar @

nzv=ones(nw,1) ;
                 @ nzv=nz(1)|..|nz(nw); nw by 1 vector
                   For each i, nz(i) is # of elements in zi(t) i=1,..,nw.
                   Let L=sumc(nzv), then L is # of orthogonal condition.@

rwv=ones(nw,1) ;
                  @ rwv=rw(1)|..|rw(nw); nw by 1 vector
                    rw(i) is defined so that
                    wi(t) is in I(t+rw(i)) i=1,..,nw @

@---------------------------------@
     @ Section B.  User can leave the variables in this section as they are,
                   and go to STEP 3.

                   Only advanced users should modify this section.@

@---------------------------------@

gradname=&GRAD2; @ Specify the name of proc that calculates
                    the gradient dgT(b)/db. @

const=1/sqrt(tend);    @ Scaling Multiplier when W0=eye(L) @
const2=1/sqrt(tend);   @ Scaing Multiplier when W0 is not eye(L) @
                       @ These scaling multipliers should be set so that
                         the value of function (vof) in nonlinear search of
                         MINQUAD.SET is near 1.  If vof is too close to 0,
                         the search will not work properly.

                         const2=1/sqrt(tend) will generally be good. @
w0flag=0;    @ scalar;
                If w0flag=0, W0=I is used as the initial weighing matrix W0.

                If w0flag=1, initial bgm is used to calculate initial W0.

                If w0flag=2, W0 in the memory is used as initial W0.

                If w0flag=3, W0 and bgm in the memory are used to give
                       the first GMM result.                       @

maxitegm=2; @ Sets maximum # of iteration over weighting matrix, W0.  Set
                 w0flag=0 and maxitegm=2 to execute usual 2-Stage GMM. @

zero=1E-2;    @ Iteration over W0 continues until the maximum
                 difference of the current and the previous W0 in
                 absolute value becomes less than 'zero', or
                 the # of iteration exceeds maxitegm. @

calwflag=4;  @ This variable is used to choose the method to calculate
               the distance matrics, W0.

               If calwflag=0, Durbin's method will be used when W0
                   is singular.

               If calwflag=1, Newey and West Method will be used when
                   W0 is singular.

               If calwflag=2, Parzen's lag window will be used when
                  W0 is singular.

               If calwflag=3, Durbin's method will be used.

               If calwflag=4, Newey and West Method will be used.

               If calwflag=5, Parzen's lag window will be used.

           Durbin's method imposes zero restrictions while Newey and West
           method dose not.  @

ordard=floor(sqrt(8190/sumc(nzv)^2)); @ Order of AR representation for
                 Durbin's method when W0 is singular @

lend=5;   @ Order of lags used for Newey and West method @

    @ See MINQUAD.SET for the following globals @
hflag=2;
dfpflag=1;
sstol=1e-25;
    @ See MAXMUM.DOC on MODULE9 of GAUSS for the following globals @
gradtol=1e-5;
btol=1e-5;
typf=1;
typb=1;

@ *************************************************************** @
@ ----------- Step 3: READING IN DATA --------------------------- @
@ *************************************************************** @

load x[115,7]=hoarding.dat ;
c=x[.,1];
dk=x[.,2];
y=x[.,4];
n=x[.,5]./1369 ;
k=x[.,6];
hc=hpfilter(ln(c),1600) ;
hi=hpfilter(ln(dk),1600) ;
hy=hpfilter(ln(y),1600) ;
hn=hpfilter(ln(n),1600) ;
hapl=hy-hn ;
clear x;

@ ******************************************************************** @
/* -------------------------------------------------------------------
          Step 4: Define the proc hu(b) that returns tend by L matrix
                             |[z(1)w(1)]'      |
                             |    |            |
                             |[z(tend)w(tend)]'|

                     where [z(t)w(t)]'=[w1(t)z1(t)',...,wnw(t)znw(t)']
----------------------------------------------------------------------- */
@ ********************************************************************** @

@ The proc must have the name "hu", and should take only one argument. @

proc hu(b);
 local dep,t,t0,w8l,beta,a,w1,w2,w3,w4,w5,w6,w7,w8,w9,w10,w11,w12,w13 ;

 dep=1-(k[2:115,1]-dk[1:114,1])./k[1:114,1] ;
 t=seqa(1,1,113) ;
 t0=0|t|114 ;
 beta=1.03^(-.25) ;
 a=y[1:115,1]./(k[1:115]^(1-b[8,1]).*(n[1:115,1].*exp(b[6,1]*t0))^b[8,1]);

 w1=b[1,1]-b[8,1]*y[2:114,1]./(c[2:114,1].*n[2:114,1]) ;
 w2=ln(a[2:114,1])-b[2,1]-b[3,1]*ln(a[1:113,1]) ;
 w3=w2.*ln(a[1:113,1]) ;
 w4=w2^2-b[4,1]^2 ;
 w5=ln(y[2:114,1])-b[5,1]-b[6,1].*t ;
 w6=w5.*t/tend ;
 w7=dep[2:114,1]-b[7,1] ;
 w8=1-beta.*(c[2:114,1]./c[3:115,1]).*((1-b[8,1])*y[3:115,1]./k[3:115,1]
     +(1-b[7,1])) ;
 w9=hy[2:114,1]^2-b[9,1]^2 ;
 w10=hc[2:114,1]^2-b[10,1]^2*hy[2:114,1]^2 ;
 w11=hi[2:114,1]^2-b[11,1]^2*hy[2:114,1]^2 ;
 w12=hn[2:114,1]^2-b[12,1]^2*hy[2:114,1]^2 ;
 w13=hn[2:114,1]^2-b[13,1]^2*hapl[2:114,1]^2 ;

retp(w1~w2~w3~w4~w5~w6~w7~w8~w9~w10~w11~w12~w13) ;
endp;

@=========== The User does not have to change the code below. ============== @

@--------------------- Print results of estimation---------------------@
proc (0)=prntrslt(x,s);

? "
--------------- Minimization Results ----------------------
Step size: " s "
Value of the objective function: " vof "
Minimiser is"  x' "
-----------------------------------------------------------";
endp;

@ ------------------------------------------------------------------------- @
@ ************************************************************************ @
@ THE PROGRAM STARTS @

chi=gmmq(gradname,tend,nzv,rwv,nw,rows(bgm),zero,maxitegm);
load bgm ;
dummy=tests(bgm) ;

output off;
@ ------------------------ END OF PROGRAM --------------------------------- @

end ;

.