*
* Example 7.7 on pp 142-143.
* Technically, this jumps a bit ahead in terms of the estimation methodology, as
* it uses GMM.
*
open data tablef5-1[1].txt
calendar(q) 1950
data(format=prn,org=columns) 1950:1 2000:4 year qtr realgdp realcons realinvs $
   realgovt realdpi cpi_u m1 tbilrate unemp pop infl realint
*
set logmp = log(m1/cpi_u)
set logy  = log(realgdp)
*
compute start=1951:3
*
* Do the regression over the full period
*
instruments constant tbilrate{0 1} logy{1 2}
linreg(inst,optimal,lags=4,lwindow=bartlett) logmp start 2000:4 resids
# constant logy tbilrate
*
* Get the full sample weight matrix, and the number of observations
*
compute wfull=%wmatrix,baset=float(%nobs)
*
* Compute the full sample X'Z matrix (this is -derivative of moment conditions
* wrt the coefficients). This defines X'Z (X = regressors, Z = instruments)
* because the regressor list is used as the first list, and the instruments as
* the second.
*
cmom(zxmatrix=xz,lastreg,nodepvar,instruments) start 2000:4
*
* Compute the S**-1 * (Z'X) x full sample covmat x X'Z S**-1 needed for the LM
* tests. (S**-1 is the full sample weight matrix).
*
compute lmv=%mqform(%xx,xz*%wmatrix)
*
* Set the change point range to roughly (.15,.85)
*
compute cstart=1957:3,cend=1993:3
*
* Set up the series for the Wald, LR and LM statistics
*
set walds cstart cend = 0
set lrs   cstart cend = 0
set lms   cstart cend = 0
do splits=cstart,cend
*
*    Figure out what the pi value is for this change point
*
   compute pi=(splits-1951:3)/baset
*
*    Get the moment conditions for the first subsample
*
   cmom(zxmatrix=m1pi,instruments) start splits
   # resids
*
*    Get the moment conditions for the second subsample
*
   cmom(zxmatrix=m2pi,instruments) splits+1 2000:4
   # resids
*
*    Do GMM on the first subsample, using (1/pi)*full sample weight matrix
*
   linreg(noprint,inst,wmatrix=wfull*(1.0/pi)) logmp start splits
   # constant logy tbilrate
   compute vpre=%xx,betapre=%beta,uzwzupre=%uzwzu
*
*    Do GMM on the second subsample, using 1/(1-pi) * full sample weight matrix
*
   linreg(noprint,inst,wmatrix=wfull*(1.0/(1-pi))) logmp splits+1 2000:4
   # constant logy tbilrate
   compute vpost=%xx,betapost=%beta,uzwzupost=%uzwzu

   compute walds(splits)=%qform(inv(vpre+vpost),betapre-betapost)
   compute lrs(splits)=%qform(wfull*(1.0/pi),m1pi)+%qform(wfull*(1.0/(1-pi)),m2pi)-$
      (uzwzupre+uzwzupost)
   compute lms(splits)=1.0/(pi*(1-pi))*%qform(lmv,m1pi)
end do splits
*
*  Graph the test statistics
*
graph(key=upleft,footer="Figure 7.7 Structural Change Test Statistics") 3
# walds
# lrs
# lms
*
* Figure out the grand test statistics
*
disp "LM" @10 *.## %maxvalue(lms) %avg(lms) log(%avg(%exp(.5*lms)))
disp "Wald" @10 *.## %maxvalue(walds) %avg(walds) log(%avg(%exp(.5*walds)))
disp "LR" @10 *.## %maxvalue(lrs)
*
* This does a sample split GMM at a specific date (here after 1980:1), using the
* theretically optimal weight matrix from Andrews.
*
nonlin b10 b11 b12 b20 b21 b22
instruments constant tbilrate{0 1} logy{1 2}
frml period1 = %if(t<=1980:1,logmp-b10-b11*logy-b12*tbilrate,0.0)
frml period2 = %if(t>1980:1,logmp-b20-b21*logy-b22*tbilrate,0.0)
compute pi=(1980:1-start)/baset
compute [symm] sw=%kroneker(%diag(||1.0/pi,1.0/(1-pi)||),wfull)
nlsystem(instruments,sw=sw,vcv) start 2000:4 period1 period2