@=================================================================================@
@ This is a GAUSS code for estimating the C-CAPM of Hansen and Singleton (1982),  @
@ using the two-step GMM, iterated GMM and continuously-updated GMM estimators.   @
@ The data is constructed from the data files in Hansen/Heaton/Ogaki GMM package. @
@ The optimization is implemented using the quasi-Newton procedure in GAUSS.      @ 
@=================================================================================@

new;
output file=gmm.out reset;
_kernel=2; @ 1 for Quadratic spectral kernel, 2 for Parzen kernel @
_const=0;
__output=0;

@ for description of the data file, see the footnote in Assignment 3 @
load d[172,9]=habit.dat;
ct_1=d[.,1]; ct=d[.,2]; ct1=d[.,3]; ct2=d[.,4];
mret=d[.,5];
y=(ct1./ct)~mret;
z=ones(172,1)~d[.,6:9]; @ vector of instruments @
nobs=rows(y); k=cols(z);

theta0={1,2}; @ initial values of parameters @
Wn=eye(k); @ initial weighting matrix=identity matrix @
j=1; do while j<=5;
 {theta,fmin,g,retcode}=qnewton(&crit,theta0);
 theta0=theta;
 ggg=hac(mom(theta),_kernel,_const)/nobs; @ HAC VCM of moment conditions @
 Wn=invpd(ggg); @ updated weighting matrix @
 M=gradp(&gen,theta); @ matrix of first derivatives @
 OM=invpd(M'*Wn*M)/nobs; @ variance-covariance matrix of parameters @
 sterr=sqrt(diag(OM)); @ standard errors of parameters @
 Jn=nobs*crit(theta); @ J-test for over-identification @
 pval=cdfchic(Jn,k-2); @ asymptotic p-value of J-test @
 if j==2;
  ? "Two-Step GMM";
  call prres(theta,sterr,Jn,pval);
 endif;
j=j+1; endo;

? "Iterated GMM";
call prres(theta,sterr,Jn,pval);

? "Continuously-Updated GMM";
theta0=theta; @ initial values of parameters @
{theta,fmin,g,retcode}=qnewton(&crit2,theta0);
ggg=hac(mom(theta),_kernel,_const)/nobs; @ HAC VCM of moment conditions @
Wn=invpd(ggg); @ updated weighting matrix @
M=gradp(&gen,theta); @ matrix of first derivatives @
OM=invpd(M'*Wn*M)/nobs; @ variance-covariance matrix of parameters @
sterr=sqrt(diag(OM)); @ standard errors of parameters @
Jn=nobs*fmin; @ J-test for over-identification @
pval=cdfchic(Jn,k-2); @ asymptotic p-value of J-test @
call prres(theta,sterr,Jn,pval);

output off;

@ Moment conditions @
proc mom(theta);
 local gi;
   gi=(theta[1]*(y[.,1]^(-theta[2])).*y[.,2]-1).*z;
 retp(gi);
endp;

proc gen(theta);
 local gn;
   gn=meanc(mom(theta));
 retp(gn);
endp;

@ Criterion function @
proc crit(theta);
 local Qn;
   Qn=gen(theta)'*Wn*gen(theta);
 retp(Qn);
endp;

proc wn_cu(theta);
local v,wn;
 v=hac(mom(theta),_kernel,_const)/nobs;
 wn=inv(v);
retp(wn);
endp;

proc crit2(theta);
 local Qn;
   Qn=gen(theta)'*wn_cu(theta)*gen(theta);
 retp(Qn);
endp;

@ Non-parametric HAC variance-covariance estimation @
proc hac(e,_kernel,_const);
local t,k,a2,bandw,j,xx,delta,kern,sig,ac,e1,e2,se,ee,ae,w,hac_vcm;
 k=cols(e);
 t=rows(e);
 ae=zeros(k,1); ee=zeros(t-1,k); se=zeros(k,1);
 j=1; do while j<=k;
  e1=e[2:t,j];
  e2=e[1:t-1,j];
  ae[j]=inv(e2'e2)*(e1'e2);
  ee[.,j]=e1-e2*ae[j];
  se[j]=meanc(ee[.,j].^2);
 j=j+1;
 endo;
 if _const==0;
  w=ones(k,1);
 elseif _const==1; @ the weight on the constant term is 0 @
  w=0|ones(k-1,1);
 endif;
 a2=4*sumc(w.*(ae.*se./((1-ae).^4)).^2)/sumc(w.*(se./((1-ae).^2)).^2);
 if _kernel==1; @  Bandwidth for Quadratic spectral kernel @
  bandw=1.3221*((a2*t)^.2);
 elseif _kernel==2; @  Bandwidth for Parzen kernel @
  bandw=2.6614*((a2*t)^.2);
 endif;
 xx=seqa(1,1,t-1)/bandw;
 if _kernel==1; @ Quadratic spectral kernel @
  delta=xx*1.2*pi;
  kern=3*((sin(delta)./delta-cos(delta))./(delta.^2));
 elseif _kernel==2; @ Parzen kernel @
  kern=(1-6*(xx.^2)+6*(xx.^3)).*(xx.<=0.5).*(xx.>=0)+
	(2*((1-xx).^3)).*(xx .<=1).*(xx.>0.5);
 endif;
 sig=moment(e,0);
 ac=zeros(k,k);
 j=1; do while j<=t-1;
   ac=ac+(e[1:t-j,.]'e[1+j:t,.])*kern[j];
 j=j+1; endo;
 hac_vcm=sig+ac+ac'; @ HAC variance-covariance matrix estimate @
retp(hac_vcm);
endp;

proc(0)=prres(b,se,jt,pv);
 print "         estimates            s.e.";
 b~se;
 print "value of objective function";;fmin;
 if k<=2; else;
 print "J-test for overidentifying restrictions";;jt;
 print "p-value of the J-test";;pv;
 endif; print;
retp;
endp;